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ABSTRACT 


A  computerized  technique  combining  backward  wave  propagation  and  an 
automatic  edge  detection  scheme  is  developed  and  tested.   The  class  of 
objects  considered  is  limited  to  those  with  edge  boundaries  since  it  is 
shown  that  a  universal  automatic  reconstruction  cannot  be  obtained  for 
all  possible  objects.   Using  samples  of  the  acoustic  diffraction  pattern 
as  input  this  technique  enables  the  computer  to  predict  the  most  likely 
locations  of  objects  and  to  produce  graphical  output  of  the  objects. 
A  discussion  is  given  on  implementing  the  diffraction  equations  on  the 
computer.   A  simplified  edge  detection  scheme  conserving  both  memory 
space  and  computer  time  is  described.   Test  results  are  presented  for 
both  computer  generated  diffraction  patterns  and  one  set  of  experimental 
data. 
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I.   MOTIVATION 

Since  Dennis  Gabor  {Ref.  1]  discovered  holography  -  recording  the 
amplitude  and  phase  of  the  diffraction  pattern  of  an  object  with  the 
help  of  a  reference  wave  and  the  subsequent  reconstruction  of  the 
original  object  from  this  hologram  -  vast  advances  have  been  made  in 
this  field.   Originally  confined  to  optics  holography  has  been  extended 
to  the  fields  of  acoustics  and  microwaves;  comparisons  between  the 
results  obtained  by  the  pioneers  in  holography  and  those  of  recent 
vintage  are  striking  with  the  promise  of  even  better  results  to  come. 
(For  a  description  of  the  development  of  holography  see  Gabor 's  lecture 
on  the  occasion  of  receiving  the  1971  Nobel  Prize  for  Physics  in 
Stockholm,  Dec.  13,  1971,  reprinted  in  [Ref.  2]). 

The  interest  in  acoustical  holography  is  motivated  by  the  search 
for  a  means  of  probing  media  which  cannot  be  investigated  by  any  other 
method.   For  example,  electro-magnetic  waves  have  a  very  short  propaga- 
tion distance  under  water  (those  that  can  be  propagated  for  a  longer 
distance  at  very  low  frequency,  have  much  too  low  a  resolution  to  be 
of  any  significance  for  imaging).   In  murky  water  optical  wavelengths 
cannot  be  used  at  all.   In  another  application  it  has  been  found  that 
X-rays  which  give  a  clear  indication  of  the  bony  structure  of  a  body 
are  not  very  good  at  making  distinction  between  two  different  kinds  of 
soft  tissue  (e.g.,  in  obstetrics).   Acoustic  waves  however  can  show  a 
difference  in  absorption  for  different  soft  tissues  and  do  not  result 
in  damage  to  the  insonified  body  at  reasonable  intensities.   Investiga- 
tion by  insonif ication  and  interpretation  of  the  reflected  sound  has 
long  been  used  in  probing  the  earth's  crust  and  in  nondestructive  testing. 
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The  application  of  holographic  concepts  brings  to  these  fields  the 
possibility  of  obtaining  information  that  not  only  indicates  the  existence 
of  an  object  of  interest  but  also  shows  many  of  its  features.   The  same 
statement  holds  for  the  advantages  of  acoustical  holography  versus  the 
well  known  echo-ranging  concepts  as  expressed  in  sonar.   The  result  of 
the  holographic  approach  is  an  image  of  the  object,  whereas  without 
holography  only  detection  of  the  object  is  realized.   This  is  not  meant 
to  imply  that  the  holographic  approach  is  the  only  one  feasible: 
Imaging  systems  incorporating  acoustic  lenses  have  also  been  developed 
with  good  results.   (For  the  development  and  state  of  the  art  of 
acoustical  holography  the  reader  is  referred  to  [Ref.  3]). 

In  general  illumination  of  a  recorded  hologram  with  a  duplicate  of 
the  reference  wave  results  in  two  reconstructed  object  wavefronts,  one 
being  an  exact  duplicate  of  the  original  object  wave  (producing  the 
"true"  image),  while  the  other  wavefront  is  the  complex  conjugate  of 
the  object  wave  (the  "conjugate"  image).   If  the  two  obtained  wavefronts 
are  not  spatially  separated  from  each  other  by  some  suitable  arrangement, 
they  can  interfere  with  each  other  which  results  in  distortion  of  the 
images  (the  so-called  twin-image  problem). 

Unfortunately  it  is  not  possible  to  realize  in  acoustical  holography 
that  most  striking  feature  of  optical  holograms,  namely  the  three- 
dimensional  sensation  experienced  by  observing  the  true  image  with  both 
eyes.   In  order  to  be  able  to  retain  the  depth  information,  the  wave- 
lengths of  constructing  and  reconstructing  beam  have  to  be  equal  or 
nearly  so.   A  hologram  can  be  considered  to  be  a  sum  of  zone  plates 
into  which  the  object  information  is  incorporated.   (See  [Ref.  4]).   It 
then  has  an  equivalent  focal  length  f  which  is  given  by 


/I  1-1 

f  -  (  v  -  u  ) 
where  v  and  u  are  object  and  reference  source  distance  respectively, 
measured  from  the  hologram  plane  along  the  line  joining  object  and 
reference  source.  See  Fig.  1-1. 

Then  if  the  linear  dimension  of  the  hologram  is  changed  by  a  factor 
m  from  X  to  X1,  (i.e.,  X1=  mX) ,  and  the  wave  length  X1  used  in  recon- 
struction differs  from  the  original  wavelength  X   by  y  =  X  A,  the 
effective  equivalent  focal  length  of  the  hologram  is  given  by 

«]  ■  ±  a-^"1 

where  f1  =  yf/m2 

v  and  u  are  the  distance  of  reconstructed  image  and  illuminating 
source  respectively,  measured  as  before,  and  the  +  sign  gives  the 
position  of  the  true  image,  the  -  sign  gives  the  position  of  the 
conjugate  image.   This  gives  rise  to  axial  magnification 

M  =    (y/m2)  (1/v2) 

A   {+  (y/m2)  [(1/v)  -  (1/u)]  +  1/u1}2 

and  transverse  magnification 

These  magnifications  are  different  from  one  another  unless  m  =  1/y 
and  thus  give  rise  to  distortion  of  the  reconstructed  image,  if  this 
condition  is  not  met  (as  is  usually  the  case).   In  acoustical  holography 
the  reconstruction  wavelength  is  optical.   The  sound  wavelengths  used 

in  the  recording  process  are  much  longer  and  y  then  is  of  the  order  of 

_3 

10   .   To  avoid  distortion  this  would  mean  reducing  the  size  of  an 

acoustical  hologram  by  this  same  factor,  i.e.  ,  a  hologram  of  lm 
diameter  would  be  reduced  to  1mm  diameter.   Naturally  one  cannot 
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Figure  1-1 
Geometry  of  General  Hologram 
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expect  to  experience  any  three-dimensional  illusion  in  looking  at  a 
hologram  of  this  small  size.   Optical  magnification  of  this  small 
hologram  introduces  longitudinal  distortion  exactly  analogous  to  that 
being  eliminated  by  the  scaling  due  to  depth  of  field  problems.   This 
longitudinal  distortion  is  one  of  the  primary  disadvantages  of  acousti- 
cal holography. 

Another  disadvantage  encountered  in  acoustical  holography  is  the 
absence  to  date,  inspite  of  intensive  research,  of  a  recording  medium 
which  gives  the  analogous  fine  grain  structure  as  that  encountered  in 
optical  holography  by  using  spectroscopic  plates.   This  shortcoming 
is  in  part  balanced  by  the  possibility  of  recording  both  amplitude 
and  phase  of  a  sound  wave  impinging  on  a  linear  detector  directly. 
Because  of  the  presence  of  linear  acoustic  detectors  (one  of  the 
fundamental  differences  of  acoustical  holography  versus  optical 
holography) ,  it  is  not  necessary  to  use  an  acoustic  reference  beam  and 
to  work  with  intensities  as  in  optical  holography.   It  might  be  argued 
that  a  direct  recording  of  amplitude  and  phase  of  a  diffraction  pattern 
is  not  really  a  hologram  (according  to  definition)  because  of  the  lack 
of  the  reference  beam. .  But  as  a  hologram  is  after  all  an  ingenious 
method  of  recording  a  complex  wavefront  using  a  medium  that  is  insensi- 
tive to  phase  and  as  the  recorded  diffraction  pattern  contains  the  same 
information  as  a  hologram,  only  in  simpler  form,  there  is  no  difference 
in  the  information  content  of  the  two  types  of  recording.   Working  with 
amplitude  and  phase  can  also  eliminate  the  twin-image  problem  referred 
to  earlier. 

Assuming  then  the  existence  of  a  record  of  the  amplitude  and  phase 

of  the  sound  diffraction  pattern  due  to  some  source  (self-emitting  or 
reflecting) ,  in  sampled  form,  whether  obtained  by  means  of  an  array  of 
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detectors  or  by  scanning  the  sound  field  or  by  any  other  means,  it  is 
possible  to  circumvent  the  difficulties  of  optical  reconstruction  (not 
only  depth  distortion,  but  also  those  due  to  degradation  because  of 
nonlinearities  of  the  film  or  cathode  ray  tube,  etc.  as  in  [Ref.  5])  by 
using  a  digital  computer  technique  to  restore  the  original  object. 
However  there  will  be  errors  introduced  by  quantization,  sampling  of 
the  complex  amplitude  and  nonunif ormities  in  the  scan  if  such  is  used. 
These  effects  are  analyzed  in  [Ref.  5].   Of  course  the  digital  computer 
is  limited  in  the  amount  of  data  it  can  handle,  whereas  the  optical 
hologram  is  comparatively  unlimited  in  the  amount  of  information  stored. 
Thus  in  increasing  the  number  of  independent  samples  of  the  wavefront 
eventually  a  point  is  reached  in  computer  reconstruction  beyond  which 
the  amount  of  storage  and  computation  time  becomes  unacceptable.   Another 
disadvantage  of  digital  reconstruction  is  that  perfect  reconstruction 
can  only  be  achieved  for  planar  objects.   The  results  of  trying  to 
reconstruct  a  three-dimensional  acoustical  object  are  identical  to 
observing  the  real  image  in  a  reconstruction  of  an  optical  hologram. 
The  shape  of  a  three-dimensional  object  can  only  be  reconstructed  in  a 
series  of  sections  through  its  real  image. 

One  of  the  greatest  drawbacks  of  digital  reconstruction  of  objects 
from  their  diffraction  pattern  is  that  the  distance  from  hologram  to 
its  generating  object  has  to  be  known.   This  problem  never  occurs  in 
the  case  of  visual  reconstruction  of  an  object  from  its  virtual  image 
because  the  focus  is  found  automatically  by  an  eye  -  brain  interaction. 
It  is  the  same  problem  as  that  experienced  in  attempting  to  bring  a 
real  image  to  a  focus  on  a  screen,  with  the  difference  that  the  move- 
ment of  the  screen  can  be  continuous  and  information  on  the  validity 
of  the  placement  is  fed  back  by  looking  at  the  result  and  interpreting 

13 


it  to  be  more  or  less  correct  than  before.   The  final  screen  location 
is  determined  by  the  compatibility  of  the  retinal  image  with  some 
object  hypothesis  stored  in  the  brain's  memory.   In  the  field  of  eye- 
brain  interaction  great  amounts  of  research  effort  have  been  and  still 
are  expended;  still  the  processes  are  not  completely  understood  to 
date.   It  is  therefore  hardly  promising  to  try  and  simulate  this  inter- 
action on  the  digital  computer.   Nevertheless  it  would  be  of  great  help 
to  find  a  method  which  brings  about  automatic  focusing  -  or  at  least 
narrows  down  that  region  in  space  in  which  the  object  might  be  present. 
Problems  can  be  envisioned  in  which  acoustical  holograms  from  quite  a 
large  number  of  different  objects  might  be  taken.   Then  one  would  not 
want  to  occupy  a  large  number  of  specialists  in  the  given  field  by 
having  them  look  at  the  results  of  each  backward  propagation  step  but 
would  rather  have  the  computer  make  a  preliminary  decision  on  the  object 
location  and  present  this  together  with  the  computed  sound  field  at 
that  location.   A  single  observer  could  then  -  using  his  pattern 
recognition  ability  -  verify  the  computer  decision  or  discard  it. 
Naturally,  computer  time  being  expensive,  the  economics  of  either 
scheme  have  to  be  considered.   In  any  case  the  computer  program  should 
take  as  little  time  and/or  memory  capability  as  possible,  an  obvious 
requirement  for  any  efficient  algorithm. 

This  then  is  the  problem  to  be  solved  in  this  investigation:   Given 
a  sampled  sound  diffraction  pattern  it  was  required  to  find  a  method 
for  a  digital  computer  to  make  a  tentative  decision  on  the  object 
location.   In  all  that  follows  a  diffraction  pattern  will  be  assumed 
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(except  where  stated  otherwise) ,  that  results  from  propagation  through 
a  medium  that  is  linear  and  not  disturbed  by  turbulence  or  convection. 
In  the  following  section  the  necessary  theory  for  this  investigation 
is  introduced. 
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II.   SCALAR  THEORY  OF  DIFFRACTION 

Diffraction  is  generally  treated  with  a  scalar  wave  theory.   This 
is  quite  applicable  for  acoustical  disturbances,  more  so  than  for 
electromagnetic  waves  which  should  rigorously  be  considered  in  a 
vectorial  theory  of  diffraction  as  their  electric  and  magnetic  vector 
components  are  not  independent  but  coupled  through  Maxwell's  equations. 

A.   DIFFRACTION  BY  A  PLANE  SCREEN 

The  solution  of  the  problem  of  diffraction  by  a  plane  screen  is 
treated  in  numerous  books,  e.g.,  [Ref.  6].   As  it  in  itself  is  not  a 
suject  of  this  investigation  only  those  results  which  are  essential  for 
the  further  treatment  are  stated  below.   Figure  2-1  shows  the  geometry 
of  the  problem. 

The  resulting  complex  amplitude  tJ  (PQ)  at  an  observation  point  PQ 
behind  an  impenetrable  screen  S,  with  an  aperture  A,  which  is  insonified 
by  a  disturbance  U  can  be  written  as 

£  (*o>  =  T-  (^'G-U^)  ds  (2-1) 

sl 

where  G_  is  the  so-called  Green's  Function,  which  is  chosen  here  to  be 
a  spherical  wave  expanding  about  the  point  PQ  and  3/3n  is  the  derivative 
in  the  direction  of  the  outward  normal  to  the  plane  of  S. 

The  difficulty  encountered  now  is  that  the  boundary  conditions 
which  should  be  applied  to  the  solution  of  Eq.  (2-1)  are  never  known 
exactly,  for  even  though  the  screen  is  assumed  to  be  impenetrable  to 
sound  (except  at  its  opening  A),  diffraction  is  known  to  take  place, 
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Figure  2-1 
Geometry  of  Diffraction  by  a  Plane  Screen 
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the  exact  expression  for  which  is  to  be  found  from  Eq.  (2-1).   To 
simplify  the  problem,  Kirchhoff  set  up  the  following  conditions  on  S^ : 

1.  Across  the  aperture  A  the  field  distribution  IJ  and  3U/3n  are 
the  same  as  they  would  be  without  the  screen. 

2.  Over  that  portion  of  S^  that  lies  in  the  geometrical  shadow. 
U  and  3U/9n  are  zero. 

It  is  noted  that  these  boundary  conditions  are  only  approximations. 
Still  Kirchhoff 's  diffraction  theory  leads  to  very  good  results  which 
are  widely  used  in  practice. 

Without  going  into  the  details  (see  [Ref.  6])  Kirchhoff 's  diffraction 
equation  for  the  plane  screen  insonified  by  a  unity  amplitude  spherical 

wave  emanating  from  P2  is, 

(2-2) 

1  exp    [ik(r2i"l_  rni)] 

U    (P   )    =  —     //  W1 ^ [cosU.roi)    -   cos(n,r21)]ds 

2jA       r2i  •  r01  *■»■ 

A 
where  (n,rQ^)  indicates  the  angle  between  the  indicated  vectors. 
Another  severe  inconsistency  of  the  Kirchhoff  boundary  conditions  is 
that  by  assuming  a  value  of  zero  for  both  U  and  3U/3ri  the  problem  is 
over-specified  and  requires  the  field  to  be  zero  over  all  space  which 
clearly  is  not  the  case.   Considerable  criticism  has  therefore  been 
raised  against  Kirchhoff s  treatment.   Sommerfeld  removed  this  diffi- 
culty by  assuming  a  Green's  function  which  makes  G_  or  3G/3n  identically 
zero  over  all  of  S,,  without  invalidating  the  development  which  leads 
to  Eq.  (2-1)  (for  details  see  {Ref.  7]).   His  result  is  incorporated 
into  the  Rayleigh-Sommerfeld  diffraction  formula: 

1  exp(ikrn , ) 

U  (PQ)  -  jX  //  U  (P,)    rQ1  01    cos  (n,701)ds  (2-3) 

A 
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where  U  (P,)  is  the  field  distribution  in  the  aperture  plane.   This  can 
be  expanded  into  (2-4) 

A  (xq-Xj)  2  +(y0-y1)2+(Z0-Z1)2 

which  shows  that  diffraction  is  a  space-invariant  phenomenon,  as  the 
result  depends  on  the  difference  between  observation  point  and  the 
point  at  which  the  disturbance  is  located. 

Diffraction  also  is  a  linear  phenomenon,  because  the  response  to  a 
sum  of  inputs  is  equal  to  the  sum  of  the  individual  responses  to  each 
input,  (within  the  limitations  of  this  study).   The  linearity  and 
invariance  of  diffraction  imply  that  the  result  Uq  due  to  some  distur- 
bance U_  impinging  on  a  diffracting  aperture  can  be  written  in  form  of  a 
convolution  integral 

Uo  =  U  *  h 
where  h_  is  the  so-called  impulse  response  of  the  diffracting  system. 
It  is  also  possible  then  to  write: 

F  {U  }  =  F  {U  *h}  =  F  {U}  •  F  {h}  =  F  {£}•  H 

where  F  is  the  symbol  for  the  Fourier  transform,  H  is  the  transfer 
function  of  the  diffracting  system  and  use  has  been  made  of  the  con- 
volution theorem. 

For  a  unity  amplitude  expanding  spherical  wave  incident  on  the 
screen  Eq.  (2-3)  becomes 

U  (P  )  -  -1  rr   eXp[jk(r^+r01)]     ,«  -     ,  ,  (2-5) 

E0M  hfI      rM.r01       cos<it,r01)  ds 
A 

which  differs  from  Kirchhoff's  result  only  in  the  so-called  obliquity 

factor  cos(n,rQ1)  compared  to  1/2  [cos(n,r"01)  -  cos(n,r21)]. 


19 


The  importance  of  Eq.  (2-3)  is  that  it  allows  the  computation  of  the 
field  at  any  observation  point  if  the  sound  distribution  in  a  plane 
across  the  direction  of  propagation  is  known. 

B.   SPATIAL  FREQUENCY  APPROACH  TO  DIFFRACTION  BY  PLANE  SCREEN 

The  following  discussion  follows  closely  [Ref.  7  and  Ref.  8]. 
Propagation  and  diffraction  being  linear,  space  invariant  phenomena  it 
is  possible  to  define  impulse  responses  and  transfer  functions  for 
different  systems.   Every  wave  distribution  can  be  expanded  into  a 
finite  or  infinite  sum  of  weighted  plane  waves.   This  is  accomplished 
by  Fourier  transforming  the  complex  disturbance  by  use  of  the  formula 

CO 

A   (fx,fy)    =  JVU   (x,y)    exp[-j2Tr(fx-x+  fy-y)]    dx  dy  (2-6) 

—•CO 

where  A  (fx,f  )  is  then  the  complex  spatial  frequency  spectrum  of  U 
and  the  spatial  frequencies  fx  and  fy  are  related  to  the  direction 
cosines  by 

fx  =  a/A  ;  fy  =  6/X;  fz=  yA 

The  direction  cosines  a,g  and  y  are  defined  to  be  the  cosines  of  the 
angles  between  the  direction  of  propagation  of  the  plane  wave  and 
the  positive  coordinate  axes.   The  equation  for  a  unit-amplitude  plane 
wave  propagating  with  direction  cosines  (a,8,y)  is  simply 

B  (x,y,z)  =  exp  [jk(ax  +B  y  +Y  z)]  (2-7) 

where  the  direction  cosines  are  geometrically  constrained  to  obey  the 
equation      a2  +  B2  +  y2  =  1 
Writing  U  (x,y,z)  as  an  inverse  transform  of  its  spectrum 

CD 

U    (x.y.z)    -//A   (fx,fy)    exp[j27T(fx.x  +   f    »y)]    dfx   df  (2-8) 
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it  is  obvious  upon  comparing  Eq.  (2-8)  and  Eq.  (2-7)  that  U^  can  be 
visualized  as  a  sum  of  plane  waves  propagating  with  direction  cosines 

a  =  x-fx,  b  =  x-fy  Y  ■  a-a2fx2)-a2fy2))  l/2 

The  complex  amplitude  of  a  plane  wave  component  is  simply  A_  (fx,fy) 
dfx  df^,  evaluated  at  (fx  =  a/A,f  =  B/A). 

In  the  following  treatment  the  impulse  response  and  transfer  function 
of  a  diffracting  aperture  will  be  found.   An  equivalent  form  [Ref .  8]  of 
the  Rayleigh  -  Sommerfeld  diffraction  formula  Eq.  (2-3)  is 


3        exp(jkr01). 


£  C>  -ir//£  <V  ^  t=™i    * 


(2-9) 


=  27    ff±  QJ'KOI     dS 

A 

8   exp(jkrQ1) 
where  KqX=  J^   [    f^   ] 

The  plane  wave  expansion  for  a  diverging  spherical  wave  of  the  form 
[exp(jkrQ^) ]/rQ^  is  known  to  be  (see  [Ref.  9]) 

oo 

exp(jkr01)    j    i 

+  (z0-Z;l)y]}  dadg 
To  find  Yq  i   the  differentiation  according  to  z^  is  performed  and  yields: 

K01  =  r%     ff   exp  {jk[(x0-x1)a+(y0-y1)B+(z0-z1)Y]}  dctdB       (2-11) 
X2 

» 
Then  « 


^ =  —  "7  exP{Jk[(x0-x1)a+(y0-y1)6  (2-10) 


2  (po>   ^  ff   H  (pl)  -//  exp{jk[(x0-x1)a+(y0-y1)6  (2-12) 

OO 

+  [  (zq-z-^)^  ]  }dad6  dx-^dy-^ 


X2 
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This  can  be  brought  into  the  form 

U  (x0,y0,z0)  =  //U1(x1,y1,z1)'t(x1,y1,z1)  (2-13) 

— OO 

.   ■  a  6 

'  //  exp{jil{Cx0-x1)a+(y0-y1)e-f(z0-z1)Y]}d^Y 

oo 

by  introducing  the  aperture  function  t(x, ,y, ,z,)  which  is  unity  inside 
the  aperture  located  in  the  opaque  screen  at  z-^  and  zero  outside  A. 

Kirchhoff  boundary  conditions  are  applied  hereto  t[  but  not  to  3U/3n. 
Separating  the  exponential  and  interchanging  the  order  of  integration 
leads  to 


00      00 


a         9, 

Uo(xo>yo>zo)  =  H    M  UtCxi>yi»zi)  exP[-J27T(Txi+ryi)]  dx1dy1 

—  oo     —  oo 

2tt  2tt  (2-14) 

•exp{-jy^rzi}    exp{Jx^(ax0+3y0+Yz0)}.  dy  d| 

where  Uj.    in  Eq.    (2-14)    replaces  Ui't.      But  now  the  inner   integral   is 
just   the  Fourier   transform  or  plane  wave  decomposition  Aj-of  JL. 

00 

d         R 
At   -  //   Ut(x1,y1,z1)    exp[-j2TT(xx1-^y1)]    dx1dy1  (2-15) 


and  the  outer  integral  can  be  written  as 

oo 

Uo(xo,y0,z0)  =  //  At(y>x)  exp{j|^-  /i-a^-B2  (ZQ-  Z±)} 


♦   expfj^Caxo+Byo)}^  d* 


(2-16) 


a  6 
Eq.  (2-16)  can  be  seen  as  an  inverse  transform  which  relates  A«(— ,— ) 

"■*  A   A 

to  UQ(x0,y0,z0) ,   so  that 

A^  e,A)  exp[j2l(z0-z1)  Sl-a*-Qz\.must   be  the  spectrum  An  (f.f) 

of  the  diffracted  and  propagated  disturbance.   Therefore  a  transfer 

function  H  (f  ,f  )  can  be  defined 
—   x'  y 

An(fx,fv) 


5  (fx'fy>  "  gf^fy)  =  -P^O-Zl^^-^y2]        <2~17> 
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The  corresponding  impulse  response  is  found  directly  from  the  Rayleigh- 

Sommerfeld  integral  Eq .  (2-3). 

1         exp(jkr01)        _ 

2  (Po)  =  IT//"  (pi>  r^ cos  (n'r01}  ds  =  (2_18) 

A 

i   °°  exp(jkrni)       _ 

=  tt  //  U  (P,)  -tCP-,)  —  cos(n,r0])  ds 

JX    —   1      x      r01  ux 


//  ^V^.  (p0'Pl)  ds 


The  impulse  response  _h  (PQ.,P..)  is  given  thus  by 

h  (p0fPl)  =  _L  exp[jkr01]cos(.  (2i9) 

~  jA      r01         01 

By  definition  then  Eq.  (2-17)  must  be  the  Fourier  transform  of 
Eq.  (2-18). 
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III.   FRESNEL  DIFFRACTION  AND  DIFFRACTION 

TREATMENT  IN  THE  SPATIAL  FREQUENCY  DOMAIN 


Having  found  the  general  form  of  the  solution  to  diffraction  by 
a  plane  screen  it  is  soon  evident  that  the  Rayleigh-Sommerfeld 
diffraction  formula  is  not  very  amenable  to  closed  form  solution. 
Even  on  the  digital  computer  one  would  expect  quite  long  computation 
times.  Therefore  approximations  are  introduced,  which  will  allow 
reducing  diffraction-pattern  calculations  to  simpler  mathematical 
operations. 

A.   FRESNEL  DIFFRACTION 

Consider  the  situation  of  Fig.  3-1.   Single-frequency  sound  is 
diffracted  by  a  finite  aperture  A  in  an  infinite  plane  opaque  screen 
which  is  located  in  the  x^-y^-  plane  at  the  origin  of  the  z-  axis. 
The  observation  region  is  assumed  to  be  parallel  to  the  x, -y, -  plane 
and  has  a  coordinate  system  (xg.yg)  attached  to  it.   This  observation 
plane  is  at  a  distance  Z  from  the  diffracting  screen. 

Using  the  Rayleigh-Sommerfeld  formula  the  field  at  a  point 

(x  ,y  )  can  immediately  be  written  as 
o  o 

00 

^O(x0,y0,Z0^  =  ff      -   ^x0'y0,xl'yl^  i^(xi>yi'0)  dxidyi         (3_1) 

—CO 

where 

I      exp(jkr01)       .  _ 
h  (xO>yo,xl'yP  =  jT  r cos(n,r01)  (3-la) 

and  the  effect  of  the  finite  aperture  A  is  incorporated  into  U  it 

being  assumed  the  U^  is  zero  outside  the  aperture  region.   This  allows 

extending  the  integration  limits  to  infinity. 
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The  first  assumption  to  be  made  is  that  the  distance  z  is  very  much 
greater  than  the  maximum  linear  dimension  of  the  aperture  A  and  that  in 
the  observation  plane  only  a  small  region  around  the  z-axis  is  of 
interest.   Again  z  is  assumed  to  be  very  much  greater  than  the  extension 
of  this  region  of  interest.   The  initial  approximations  introduced  are 
called  the  paraxial  approximations  because  they  hold  only  close  to  the 
axis  of  propagation.   Using  this  assumption  the  obliquity  factor  is 
approximated  by 

cos  (n,r01)  =1 

which  is  good  to  within  5%  if  the  argument  of  the  cosine  is  not  greater 
than  18°. 

The  same  assumption  can  also  be  used  to  replace  the  radius  r~,  by 
the  distance  z  in  the  denominator  of  Eq.  (3-la)  but  not  in  the  exponent, 
since  the  resulting  error  would  be  multiplied  by  the  large  number  k, 
generating  intolerable  phase  errors.   The  impulse  response  has  now  been 
approximated  by: 

-  (x0'y0'xl'yl^  "  ^T~   exP(Jkr0l)  (3~2) 

j  Az 

Further  simplification  is  possible  by  using  the  Fresnel  approxima- 
tion.  This  consists  of  expanding  the  radius  r^,  in  the  exponent,  which 
is  given  exactly  by 

r01  ■  /(xq-Xj^)  z+(y0-y1)  z+zz 


=  z  /i+(XQ-xi)/+  (yp-yip 


in  a  binomial  expansion  and  assuming  that  rn  is  given  adequately  by 
the  first  two  terms  in  the  series.   Then 


r0,    £ZU  +  1   (X0-XT)2+  \  Ctt=&>2] 
Ui  2  Z  /  Z 
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The  impulse  response  can  then  be  written  as 

h  (x0,yo,x1>7l)  =  eX^fz)  exp  {JI|  [(x0-x1)2+(yo-yi)2]}       (3-3) 

The  region  for  which  this  approximation  is  exact  enough  is  called  the 
region  of  Fresnel  diffraction.   The  accuracy  of  this  approximation 
hinges  on  the  relative  sizes  of  the  aperture  and  of  the  observation 
region,  the  propagation  distance  z  and  the  wavelength  A.   One  might 
require  that  the  maximum  phase  change  contributed  by  the  third  order 
term  in  the  series  expansion  should  be  rather  small,  say  much  less 
than  one  radian.   This  is  the  case  if 

Z3>>  4X  ^xO-xl>2+<yO-^2^SAX  (3-4)' 

For  an  observation  region  and  diffraction  region  of  about  50X  by  50X 
this  would  require  the  distance  z  to  be  much  greater  than  about  250-300A. 
In  fact  the  requirement  of  Eq.  (3-4)  is  too  restrictive  because  all 
that  is  required  is  that  the  higher  order  terms  do  not  change  the  value 
of  the  superposition  integral  of  Eq.  (3-1)  by  an  appreciable  amount. 
For  distances  z  smaller  than  those  given  by  Eq.  (3-4)  the  oscillations 
of  the  quadratic  phase  factor  given  by  Eq.  (3-3)  will  generally  be  so 
very  rapid  as  the  integration  is  performed  that  they  virtually  cancel 
out  and  only  those  points  of  the  aperture  where  x0=x,  and  yn=y-i  will 
give  a  noticeable  contribution  to  the  integral.   Near  these  points  of 
so-called  stationary  phase  the  magnitude  of  higher  order  terms  will 
normally  be  negligible.   The  principle  of  stationary  phase  is  discussed 
in  [Ref.  6].   None  the  less,  Fresnel  diffraction  cannot  be  assumed  to 
hold  for  distances  which  for  the  given  example  might  be  less  than  about 
100X  (the  paraxial  approximation  is  implied  also! )  . 
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If  one  is  satisfied  as  to  the  validity  of  the  Fresnel  approximation 
for  a  given  case,  the  superposition  integral  can  be  expressed  in  the 
form  of  the  following  equation 

tt    r  \        exp(jkz)  jk  2        2 

U   (x0,y0)   =  — : exp    [-_  (x0z+yQz)J 

u  jAz  2z        u        u 


•Jk    /..   2,       2m    _.._    r    ^2lT 


(3-5) 

*  ./v    ut(xi»yi*  exp  ^ii  (xiz+yiz^  exP  E-Jxi-  (xoxi+yoy^ 

dx^dy^ 


by  expanding  the  quadratic  terms  in  the  exponent.   This  can  be  seen 
to  be  (aside  from  multiplicative  phase  factors  and  attenuation  indepen- 
dent of  (x,,y,))  the  Fourier  transform  of  IJ  (x,  ,y, )  •  exp[2k-(x,  2+y,2)  ]  , 
evaluated  at  the  spatial  frequencies  f  =  x~/\z   and  f  =  yr./xz,  which 
assures  correct  scaling  in  the  observation  plane.   As  such  the  integral 
Eq.  (3-5)  is  amenable  to  digital  computer  solution  using  a  Fast  Fourier 
Transform  algorithm. 

It  is  also  possible  to  transform  the  impulse  response  of  Eq.  (3-3) 
because  it  is  space-invariant  (it  depends  only  on  the  differences 
(xq-x-^)  and  (yg-yi) )  to  get  a  transfer  function: 

H  (f   f  )  =  exp(jkz)  exp[-JTrXz(f  2+f  2)],  (3-6) 

A  y  A   y 

which  relates  the  spatial  frequencies  in  the  diffracting  and  in  the 
observation  plane.   This  is  immediately  observed  to  be  an  approximation 
to  the  exact  transfer  function  of  Eq.  (2-17). 
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B.   FRAUNHOFER  DIFFRACTION 

For  more  severe  restrictions  on  the  propagation  distance,  it  is 
possible  to  simplify  the  diffraction  equation  even  more.   If  z  is  large 
enough   to  assume  that  the  quadratic  phase  factor  exp[jk(x,2+y  2)/2z]  is 
approximately  unity,  then  the  diffraction  equation  reduces  to  a  straight- 
forward Fourier  transform  (again  evaluated  at  the  scaled  values  of 
spatial  frequency)  with  an  additional  multiplication 

U  (x0,y0)  =  £XP^kz)  exp  [jk(x02+y02)/2Z] 

(3-7) 

X    *!    £   (x,,y,)    exp    [-j2tt  (x^+y^)  /\z]    dx-^y-j^ 

—CO 

This  equation  is  called  the  Fraunhofer  diffraction  formula.   Unfortu- 
nately, the  distance  z  would  have  to  be  very  large  indeed  for  Eq.  (3-7) 
to  be  a  good  approximation  to  the  actual  physical  situation.   For  the 
dimensions  used  previously  to  find  the  minimum  distance  for  Fresnel 
diffraction  z  would  have  to  be  very  much  larger  than  15000X.   Neverthe- 
less, the  Fraunhofer  diffraction  equation  can  be  very  useful  in  optics, 
because  it  fits  the  situation  at  the  focal  plane  of  a  well-corrected 
lens;  it  is  of  less  use.  in  acoustical  diffraction  problems  because  of 
the  difficulty  of  building  good  acoustical  lenses.   In  the  context  of 
this  investigation  it  can  be  used  to  check  the  results  obtained  by  use 
of  the  Fresnel  equation  because  Fraunhofer  diffraction  is  the  limit 
of  very  large  propagation  distance  z  in  the  Fresnel  diffraction  equation 
and  therefore  a  Fourier  transform  of  the  diffracting  aperture  would  be 
obtained. 
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The  transfer  function  fcq.-    (3-6))  approach  to  Fresnel  diffraction  can 
also  be  used  in  the  Fraunhofer  diffraction  region,  but  one  would  hardly 
under  normal  computational  circumstances  try  to  implement  the  Fresnel 
(or  Fraunhofer)  diffraction  equation  by  using  the  transfer  function  in 
the  spatial  frequency  domain.   This  will  take  longer  computation  time 
by  requiring  a  forward  and  an  inverse  Fourier  transform,  plus  a  matrix 
multiplication  whereas  in  programming  the  Fresnel  equation  in  the  space 
domain  one  needs  only  a  forward  Fourier  transform  with  additional 
matrix  multiplications. 

C.   SPATIAL  FREQUENCY  APPROACH 

If  it  has  been  decided  that  the  use  of  longer  computation  time  is 
feasible,  then  the  implementation  of  the  spatial  frequency  diffraction 
equation  (Eq.  (2-26))  is  a  much  more  powerful  tool.   (Indeed,  for 
large  array  sizes  the  use  of  the  spatial  frequency  approach  may  be  more 
economical,  as  the  number  of  computations  for  a  Fast  Fourier  Transform 
increases  like  N'log^,  whereas  N2  multiplications  are  required  for  each 
of  the  quadratic  phase  factors).   By  using  it  not  only  the  Fresnel 
requirement  on  propagation  distance  but  also  the  paraxial  approximations 
are  circumvented.   Any  inaccuracies  in  the  results  obtained  by  the 
spatial  frequency  approach  are  thus  due  only  to  the  use  of  Kirchhoff 
boundary  conditions  on  U  or  on  3U/an.   But  as  Kirchhoff' s  diffraction 
theory  as  used  in  the  Fresnel-Kirchhoff  diffraction  formula  or  in  the 
Rayleigh-Sommerfeld  diffraction  formulation  (the  relative  accuracy  of 
both  equations  being  still  a  subject  of  active  research)  gives  results 
which  are  in  excellent  agreement  with  practice,  any  reservation  against 
the  use  of  Kirchhoff's  boundary  conditions  can  only  be  very  slight. 
Moreover  they  are  applied  in  the  Fresnel  equation  also. 
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To  program  diffraction  by  a  plane  screen  using  the  spatial  frequency 
approach  one  first  multiplies  the  insonifying  disturbance  by  the  aperture 
function,  performs  a  forward  Fourier  transform  and  multiplies  each 
spatial  frequency  component  by  the  appropriate  phase  delay  as  given  by 
Eq.  (2-26).   Inverse  transforming  then  results  in  the  diffraction 
pattern  in  the  observation  plane. 
In  shorthand  notation: 

U  (x0,y0)  =  F-l  {H  (fx,fy)  •  F  {UxCxi,  yi) • t  (xljyi)}}  (3-8) 

where  F  and  P~*  symbolize  forward  and  inverse  Fourier  transform  and 

H  (f  ,£J  is  given  by  Eq.  (2-26) 
x.      y 

H  (fx,  fy)  =  exp  [j-^p  /l-(Xfx)z-(Afy)z]  (2-26) 

JJ  is  the  insonifying  wave  complex  amplitude  and  Jt  is  the  aperture 

function. 

The  aperture  function  is  usually  taken  to  be  unity  inside  the 
aperture  and  zero  outside  but  there  is  really  no  necessity  for  doing 
that.   The  aperture  might  also  be  characterized  by  the  introduction  e.g., 
of  a  constant  phase  shift  or  a  quadratic  phase  shift  (like  that  due  to 
a  lens)  or  a  combination  of  amplitude  attenuation  and  phase  shifts. 
Of  course  the  effect  of  the  aperture  could  also  be  handled  in  the 
spatial  frequency  domain  by  using  the  convolution  theorem  of  linear 
systems  analysis  and  writing 

F  (U1(x1,y1).t(x1,y1)}  -  Ax(fx,fy)  *  T  (fx,fy) 
where  A  (f  ,f  )  is  the  spatial  frequency  spectrum  of  U,,T  (f  ,f  )  is 

—     X    y  r  i  J         r  — j_»_  -     x>  y' 

the  spatial  frequency  spectrum  of  t^  and  *  represents  the  convolution 
operation. 
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This  explanation  may  be  helpful  in  visualizing  the  effect  of  the 
aperture  in  broadening  the  angular  spectrum  of  the  disturbance  but  would 
not  be  used  on  the  computer  if  it  can  be  avoided. 

The  transfer  function  H  (f  , f  )  has  to  be  investigated  more  closely. 


To  recapitulate 


.2ttz 


H  (fx,fy)  =  exp  [j^p  A-Uf^-iXfy)*]  (2-26) 

When   the  direction   cosines   a=X*f     and  g  =A*f      satisfy  ct2+B2<l   (case   of 

x  y 

homogenous  waves)  the  effect  of  propagation  over  a  distance  z  is  just 
a  change  in  the  relative  phases  of  the  spectral  components:   Since  each 
component  travels  at  a  different  direction  it  has  to  travel  a  different 
distance  to  reach  the  observation  plane  and  thereby  relative  phase  shifts 
are  introduced.   But  when  a     +32>1  (case  of  inhomogenous  or  evanescent 
waves)  the  square  root  in  Eq.  (2-26)  is  imaginary  and 

H  (fx,fy)  =  exp  [-2*2   /Afx)  *+  (Afy)z-1]   =  exp  J-pz] 

where  p  is  a  positive  real  number.   These  wave  components  which  travel 
perpendicular  to  the  direction  of  propagation  are  severely  attenuated 
by  even  small  propagation  distances.   For  propagation  distances  greater 
than  a  few  wavelengths  the  influence  of  these  evanescent  waves  is 
entirely  negligible  on  the  diffraction  pattern  and  it  is  possible  to 
band-limit  the  transfer  function  to 


exp[j2liz  /l-(Afx)2-(Xfv)2]  for  f  2+f2<  — 
H  (fx,fy)         *  y  y        X2 

0  elsewhere 


(3-9) 
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D.   COMPARISON  OF  FRESNEL  EQUATION  AND  SPATIAL  FREQUENCY  APPROACH 

A  comparison  between  the  two  possible  ways  of  implementing  the 
diffraction  problem  is  now  possible. 

The  time  needed  for  one  propagation  step  using  the  Fresnel  equation 
(one  Fourier  transform,  two  multiplication  operations  on  a  two-dimen- 
sional matrix  -  one  before,  one  after  the  Fourier  transform)  will  be 
less  (for  moderate  array  size)  than  that  consumed  by  one  propagation 
step  using  the  spatial  frequency  approach  (forward  Fourier  transform, 
matrix  multiplication,  inverse  transform).   The  Fresnel  equation 
cannot  be  used  for  small  propagation  distances  because  of  its  inherent 
approximations,  which  do  not  occur  in  the  spatial  frequency  method. 
The  spatial  increment  between  array  points  increases  with  propagation 
distance  in  the  use  of  the  Fresnel  equation  due  to  the  necessary 
evaluation  of  the  Fourier  transform  at  fx=x0/^z,  fy=y0/^z,  whereas 

it  is  constant  .when  the  spatial  frequency  approach  is  used.   This  leads 
to  the  fact  that  it  is  possible  to  use  the  Fresnel  equation  program  for 
very  large  propagation  distances  whereas  the  spatial  frequency  method 
is  limited  in  the  propagation  distances  to  be  used  because  the  diffrac- 
tion patterns  from  "ghost"  objects  (which  are  postulated  by  the  Discrete 
Fourier  Transform)  start  overlapping  with  the  useful  diffraction  pattern. 
These  effects  are  described  in  more  detail  later  in  the  section  on  the 
Discrete  Fourier  Transform.   Reduction  of  this  unwanted  and  damaging 
effect  is  possible  but  only  by  using  a  larger  array  size  which  can 
increase  the  computation  time  considerably.   The  last  aspect  under 
which  the  two  methods  will  be  compared  is  that  of  the  validity  of  the 
results  after  backward  propagation  of  the  observed  wavefront  when  the 
propagation  distance  differs  by  a  small  amount  from  the  actual  distance. 
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Suppose  a  diffraction  pattern  is  observed  at  a  distance  z  from  a  plane 

object.   Then,  on  propagating  the  wave  back  a  distance  z-dz  (an  error 

which  might  occur  in  an  automatic  reconstruction  scheme)  using  the 

transfer  function  approach  one  would  expect  to  obtain  the  same 

diffraction  pattern  as  would  appear  on  forward  propagation  of  the 

original  diffracted  wavefront  at  the  object  by  a  distance  dz.   It  is 

then  possible  (by  propagation  backward  using  a  propagation  distance  dz) 

to  obtain  the  actual  diffracting  object.   The  case  is  quite  different 

when  the  Fresnel  equation  is  used.   Here  when  the  propagation  distance 

in  reconstructing  differs  by  a  small  amount  from  the  true  one,  the 

obtained  "image"  is  only  a  very  coarse  approximation  to  the  diffraction 

pattern  at  a  distance  dz  in  front  of  the  diffracting  aperture.   The 

actual  diffracting  plane  and  that  just  reconstructed  are  related  only 

by  the  fact  that  they  give  rise  to  the  same  diffraction  pattern  at  a 

distance  z  and  z-dz  respectively  when  the  Fresnel  equation  is  used. 

In  this  second  case  it  is  not  possible  by  another  small  propagation 

step  of  the  correct  difference  dz  to  reach  the  correct  result:   The 

Fresnel  equation  does  not  hold  true  for  small  propagation  distances. 

One  would  have  to  go  back  to  the  original  diffraction  pattern  and  to 

try  again  from  there.   This  either  requires  storing  the  original 

diffraction  pattern  (because  in  the  use  of  the  Fourier  transform  the 

original  data  are  replaced  by  the  transformed  components),  which  being 

a  two-dimensional  matrix  of  complex  quantities  with  a  rather  large 

size  requires  a  considerable  amount  of  additional  core  memory  or 

demands  another  propagation  step  for  the  restoration  of  the  original 

pattern.   Thus  the  use  of  the  Fresnel  equation  in  an  automatic 

reconstruction  scheme  would  take  more  computation  time  or  core  memory 

than  the  spatial  frequency  method  and  become  uneconomical.   (If  the 
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propagation  distance  is  known  precisely  beforehand  this  reversal 
of  relative  time  consumption  between  the  two  methods  would  not  occur) . 
This  last  aspect  settled  the  question  which  algorithm  should  be  used  in 
favor  of  the  transfer  function  approach  in  that  it  is  desirable  in  a 
restoration  problem  to  reach  a  result  which  is  as  close  to  the  actual 
disturbance  at  a  particular  distance  as  possible  even  though  the 
distance  might  be  slightly  wrong.   Moreover  it  is  much  more  difficult 
to  handle  a  diffraction  problem  by  the  Fresnel  equation  if  two  diffracting 
objects  are  located  in  planes  which  are  displaced  in  depth,  due  to  the 
self-scaling  of  the  Fresnel  equation;  it  is  impossible  to  do  if  the 
distance  between  these  objects  is  less  than  that  required  for  the 
Fresnel  equation  to  hold. 

E.   BACKWARD  PROPAGATION 

How  'backward  propagation'  is  accomplished  can  now  be  described. 
The  conceptually  more  satisfying  spatial  frequency  approach  is  treated 
first  because  it  is  used  in  the  actual  program.   Consider  a  given 
diffraction  pattern  Uq  (xq^q)  at  a  distance  z  from  an  object  U-^Cx^.yi) 
as  in  Fig.  3-1.   Then  from  the  diffraction  equations, 

Uq  (x0,y0)  =  F"1  {H  (fx,fy)-F  {t^  (xx ,yi) } }  (3-8) 

To  find  U^  (x^,y^)  one  would  first  perform  a  Fourier  transform  on  Uq, 
multiply  the  result  by  H""1  and  then  inverse  transform  the  result 

S2  (x1>yi)  =  F"1  Cff-l(fx,fy).F  {UqCxq.Vq)}}  (3-10) 

where  H-*  is  defined  by 

H'H"1  "  1  (3-11) 
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The  symmetry  between  Eq.  (3-8)  and  (3-10)  is  obvious.   The  reverse 
transfer  function  H   can  easily  be  found  from  the  definition  Eq.  (3-11) 
and  Eq.  (3-9)  which  defines  the  transfer  function  H. 

H"1  =  exp  [-j2l5.  /l-(Xfx)MAf  )z  ].  (3-12) 

Observe  that  the  same  result  would  have  been  obtained  by  mechanically 
inserting  (-z)  as  the  propagation  distance  in  the  equation  for  the 
forward  transfer  function,  where  the  minus-sign  might  be  taken  to 
signify  'backward  propagation'.   Another  observation  is  in  order  here: 
Suppose  H  were  not  taken  to  be  band-limited  to  those  spatial  frequencies 
which  give  a  real  square  root.   Then  at  least  H~  must  be  band-limited; 
else  for  those  frequencies  which  give  an  imaginary  square  root 

H"1  (fx»V  =  exP  IHT  /(Xfx)2  +  (Xfy)2-1]  =e+pZfor  f^+f^-I 

and  an  exponential  growth  would  be  postulated,  which  leads  to  conceptual 
problems  (for  complete  treatment  of  this  subject  see  [Ref.  8]).   Apart 
from  these  conceptual  problems:  the  computer  will  not  handle  exponentials 
with  real  positive  arguments  which  are  too  large  (as  would  be  in  the  case 
for  non-bandlimited  H   used  in  a  propagation  over  more  than  a  minuscule 
distance).   The  frequency  components  of  the  angular  spectrum  correspond- 
ing to  evanescent  waves  in  a  diffraction  field  that  is  several  wave- 
lengths away  from  the  object  would  be  very  small  in  amplitude.   There 
is  a  good  chance  that  they  would  be  prone  to  round-off  errors  on  the 
computer  and  in  an  actual  data  acquisition  system  would  be  lost  in  the 
unavoidable  noise  anyway.   Trying  to  work  with  a  non-bandlimited 
reverse  transfer  function  then  is  bound  to  fail.   But  it  must  always 
be  kept  in  mind  that  the  evanescent  waves  are  physically  existent  and 
carry  information  about  those  object  features  which  are  smaller  than 
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a  wavelength.   Setting  them  to  zero  destroys  information  irretrievably 
and  therefore  degrades  system  performance.   Figures  3-2  through  3-5 
give  an  example  of  this  image  degradation.   Figure  3-2  shows  the 
amplitude  distribution  of  the  object,  which  is  32  sampling  points  wide 
and  16  sampling  points  deep,  centered  in  a  64  sampling  points  square 
array.   If  the  emission  at  all  sampling  points  is  assumed  to  be  in 
phase  (as  in  a  specular  object)  and  the  sampling  points  are  one  half 
wavelength  apart  the  result  after  transforming  and  reconstruction 
neglecting  evanescent  waves  is  given  in  Fig.  3-3.   The  degradation  is 
not  very  severe.   For  in-phase  emission  and  sampling  every  0.3  wave- 
lengths the  result  of  Fig.  3-4  is  obtained  which  shows  much  lower  slope 
and  higher  overshoots.   It  still  gives  a  good  indication  of  the  actual 
shape  of  the  original  object.   Figure  3-5  is  the  result  of  sampling  0.5 
wavelengths  apart  and  assuming  a  random  phase  distribution  for  the 
object  to  be  reconstructed  (as  in  a  diffuse  object).   As  a  diffusely 
emitting  object  will  have  more  energy  in  high  frequency  components 
one  would  expect  the  loss  of  information  to  be  worse  than  in  the 
specular  case;  still  the  degree  of  degradation  as  shown  in  Fig.  3-5  is 
rather  unexpected.   For  experimental  results  of  backward  propagation 
using  the  spatial  frequency  approach  see  [Ref.  10]. 

From  the  Fresnel  equation  U,  (x,  ,y,)  can  be  found  by  forming  the 
complex  conjugate  of  U^  (xQ,y0)  and  propagating  U*  a  distance  (+z) 

U2(x2,y2)  =  eX?.[lkz)    exp  fjJi(x22+y22)]  (3_13) 

to 
.  //  Uq*  (x0,y0)  exp[jJi(x024y02)]  exp  [-jf^x^+y,^) ] 

dxody0 
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Conceptually,  the  term  'backward  propagation'  is  not  applicable  to  this 
operation;  it  is  rather  a  forward  propagation  of  a  different  though 
closely  related  diffraction  pattern.   The  result  of  this  operation  is 
the  original  object  with  the  positive  sense  of  the  x,  and  y, -  axes 

reversed,  i.e.,  U?(x-,y«)  =  U, (-x, ,-y1 ) ;  this  usually  is  no  reason 

for  not  using  this  approach,  since  the  inversion  is  easily  corrected. 
The  main  reason  for  discarding  the  Fresnel  equation  approach  to  back- 
ward propagation  is  that  it  is  not  possible  to  use  repeated  small 
propagation  steps  in  the  Fresnel  equation  due  to  its  inherent 
approximations . 
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IV.   DISCRETE  AND  FAST  FOURIER  TRANSFORM 

Diffraction  pattern  computations  can  be  accomplished  in  various 
ways,  but  if  one  of  the  approaches  as  outlined  in  the  last  section  is 
chosen  (for  the  sake  of  simplicity  and  short  computation  time)  it  will 
be  necessary  to  take  the  Fourier  transform  of  some  input  array.   Since 
it  is  possible  to  work  only  with  discrete  samples  on  the  digital 
computer  a  discrete  Fourier  transform  has  to  be  used. 

A.   DISCRETE  FOURIER  TRANSFORM 

Suppose  a  function  f(x)  is  represented  on  the  interval  [0,X]  by  N 

samples  f(n«a)  where 

X 
a  -  jj      and 

0  <_     n  _<  N  -  1 

The  Discrete  Fourier  Transform  (DFT)  of  the  sequence  f(o),  f(a),*"' 

f((N-l)«a)  is  defined  to  be  the  sequence  F(o),  F(Q),-"    F((N-l)«fi) 

where 

N-l 
F(kfi)    =   Z    f (n-a)    exp(-jn-a-k'ft)  0   <k   <   N-l  (4-1) 

n=o 

and  q  =   2tt   =  AIL.  (4-la) 

X         N-a 


,   N-l 
f(n-a)    =  ±     Z   F(k-fJ)   exp(jkttna)  (4-2) 


The    Inverse   DFT   is    given  by 

1  N"1 
-     Z 

N   k=o 
Remarks : 

1.   The  similarity  between  the  DFT  as  defined  above  and  the 
exponential  form  of  the  Fourier  series  expansion  is  obvious.   There 
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f(t)  =i    Z     F(nu)  expCjn^t)       (u^) 

T  n=  _oo  1 

and  t/2 

F(nu)  =      f(t)  exp(-jnwt)  dt 
-T/2 

In  the  case  of  the  DFT  the  integral  is  replaced  by  a  sum  as  there  are 
only  discrete  samples  of  f(x). 

2.  The  DFT  gives  as  many  frequency  components  as  there  are  samples 
of  the  function  to  be  transformed. 

3.  In  order  to  make  the  definition  of  the  DFT  in  form  of  a  Fourier 
series  consistent  it  has  to  be  assumed  that  the  sampled  function  is 
periodic  with  period  X.   Otherwise  it  could  not  give  rise  to  discrete 
frequency  components.   Thus  each  object  must  be  taken  to  have  'ghost 
objects'  of  equal  shape  located  next  to  it  at  a  periodic  distance  X 

in  both  horizontal  and  vertical  directions.   Analogously,  the  sequence 
of  spectral  components  must  be  assumed  to  be  periodic  outside  the 
actually  computed  spectral  period. 

4.  Even  if  the  sampled  function  values  are  real,  the  sequence 
of  frequency  components  is  in  general  going  to  be  complex. 

5.  The  spacing  between  frequency  components  is  the  reciprocal 

of  the  period  X,  multiplied  by  2tt .   As  the  spatial  frequency  approach 
to  diffraction  was  shown  to  be  bandlimiting  in  the  last  section,  it 
is  possible  to  compute  the  minimum  sampling  rate,  which  is  necessary 
to  avoid  overlapping  of  spectral  components.   In  the  spatial  frequency 
domain  the  spectrum  of  the  transformed  object  is  repeated  with  a 
frequency  period  of  1/a.   If  a  is  chosen  to  be  large  the  individual 
spectra  will  begin  overlapping;  it  is  then  not  possible  to  extract  the 
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original  object  from  this  distorted  spectrum.   This  effect  is  called 
'aliasing'.   To  recapitulate: 


exp 


(j2£z  /l-(^fx)2.(Xfy)2)     (fx2+y<^T 


H  (fx,fy)  = 

0  elsewhere 

This  means  that  in  the  spatial  frequency  plane  H  (f  ,f  )  is 

bandlimited  to  those  frequencies  which  lie  inside  the  circle  of  radius 
1/X  centered  at  the  origin.   (See  Fig.  4-1)   If  efficient  sampling 
is  desired  (meaning  sampling  that  is  not  higher  than  the  Nyquist  rate) 
the  sampling  scheme  becomes  very  complicated  even  for  such  a  simple 
band-limited  frequency  region  as  a  circle.   This  can  be  avoided  at  the 
cost  of  slightly  higher  sampling  rate  (about  12%,  see  [Ref.  11])  by 
assuming  the  band-limited  region  to  be  the  smallest  square  that  circum- 
scribes the  actual  circle.   Using  the  Nyquist  sampling  criterion  on 
this  new  band-limited  region  leads  to  the  maximum  allowable  spacing 
between  sampling  points,  namely  a  =  X/2.   This  is  found  the  following 
way: 

For  the  fx  -direction:    f^j^x  =  J  =   Afx  N/2 

Afv  = 


N«X    2tt 

from  Eq.  (4-la)  a  .  2w  m   2«2tt 

X     N-X 

then     X  =  — ;  a  -  -  =  - 
A  similar  result  holds  analogously  for  the  f   -direction. 
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Figure  4-1 
Propagation  Transfer  Function  in  the  Spatial  Frequency 
Domain,  Showing  Band-limitation  at  Nyquist  Rate  Sampling. 


►  f. 


Figure  4-2 
Propagation  Transfer  Function  in  the  Spatial  Frequency  Domain, 
Showing  Band-limitation  at  Twice  Nyquist  Rate  Sampling. 
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Now,    suppose   for  some   reason  it  would  be   desirable   to   introduce   a 
guard  band  in  the  spatial   frequency   domain   or   to   sample   at   a  higher  than 

the  Nyquist   rate,    e.g., 

1        X.  X 

3     "  4<a=2 

without  increasing  the  total  number  of  samples. 
Then  f  MAY  =  JU  =  1   •  (See  Fig.  4-2) 

But  the  Discrete  Fourier  Transform  will  space  the  calculated  frequency 
components  evenly  from  -fxmax  to  +  f^ax,  and  from  -fymax  to  +f  max 

and  all  the  spectral  components  outside  the  bandlimiting  circle  will  be 
set  to  zero  on  propagation.   Therefore  the  (necessarily  discrete) 
propagation  algorithm  will  have  to  work  on  fewer  frequency  components 
than  when  the  Nyquist  rate  is  used.   The  results  will  be  in  error 
because  on  inverse  transforming  fewer  spatial  frequency  components 
are  available.   There  is  then  a  trade-off  between  fine  sampling  in 
space  and  a  coarsening  of  results  of  the  propagation  algorithm  which 
has  to  be  taken  into  account. 

Of  course  if  the  object  spectrum  itself  is  band-limited  to  a  band 
which  is  narrower  than  that  introduced  by  the  propagation  transfer 
function,  the  sampling  rate  could  be  adjusted  accordingly  to  a  coarser 
sampling  (see  [Ref.  12]).   Since  the  object  is  assumed  to  be  unknown, 
this  possibility  cannot  be  used  here. 

The  definition  of  the  Discrete  Fourier  Transform  is  readily  extended 

to  two  (or  more)  dimensions.   In  particular  for  two  dimensions: 

N-l  M=l 

F  (kft.lA)  =   Z    Z  f(na,mb)  exp  [-1  (nakfl  4-  mblA)] 
n=o  m=o  r   J 

for  0  <  k  <  N-l  and  0  <_  I  <_  M-l  (4-3) 

where  0  =  jj*  and  A  =  ^;  a=  |  ;  b  =  £ 
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Then 

N-l  M-l 

f(na,mb)  =  W7m     e   z  F(kO,lA)  exp[ j (nakfi+mblA) ]  (4-A) 

iN  ra  k=0  1=0 

The  postulated  periodicity  of  the  object  to  be  transformed  is  very 
unfortunate  for  a  system  that  is  to  simulate  imaging  and  diffraction. 
Consider  the  one-dimensional  case  as  shown  in  Fig.  4-3.   An  object  and 
its  two  nearest  'ghost  objects'  are  shown.   It  is  well  known  that 
diffraction  means  a  spreading  out  of  the  diffracted  wavefront.   The 
waveform  of  Fig.  4-3  diffracted  and  propagated  out  some  distance  might 
look  like  Fig.  4-4.   Now  at  an  even  larger  propagation  distance  the 
spreading  out  might  be  so  severe  that  the  wavefronts  due  to  the  object 
of  interest  and  those  due  to  the  postulated  ghost  objects  to  either  side 
of  it  will  begin  overlapping  and  interfering  constructively  or  destruc- 
tively as  in  Fig.  4-5.   Naturally  the  waveform  of  Fig.  4-5  looks  quite 
different  towards  the  boundaries  of  the  sampling  period  from  the  one 
that  would  be  obtained  without  any  ghost  objects  and  in  fact  is  obtained 
and  observed  in  a  physically  existing  system  as  in  Fig.  4-6.   The  observed 
waveform  of  Fig.  4-6  sampled  and  fed  to  the  computer  to  be  propagated 
back  the  correct  distance  using  a  DFT  routine  cannot  reproduce  exactly 
the  central  of  Fig.  4-3,  because  of  the  limitations  of  the  DFT  routine. 
It  is  this  fact  that  gives  rise  to  the  limitation  in  propagation  distance 
of  the  transfer  function  approach  to  diffraction  that  was  mentioned  in 
the  last  section.   There  is  an  easy  wasy  to  overcome  this  limit  by 
Increasing  the  guard  space  between  the  desired  object  and  its  ghost 
images  by  increasing  the  array  size,  if  the  implied  increase  in  core 
memory  needed  and  associated  computing  time  can  be  afforded.   In  this 
investigation  it  was  decided  that  a  64  by  64  element  matrix  was  as 
large  as  could  be  afforded  without  running  into  time  and  memory  problems. 
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On  this  basis  the  largest  useful  propagation  distance  was  found  to  be 
about  200  wavelengths  for  an  object  aperture  with  a  linear  dimension  of 
the  order  of  15A. 

In  image  processing  it  is  often  desirable  to  have  both  the  center 
of  the  input  object  and  the  zero  frequency  component  (or  undiffracted 
or  d.c.  component)  of  the  spatial  frequency  spectrum  centered  in  the 
array.   When  the  DFT  relationship  of  Eqs.  (4-3)  and  (4-4)  is  used  the 
zero  frequency  component  appears  in  one  corner  of  the  transformed  array. 
Andrews  [Ref.  13]  recommends  multiplying  the  array  to  be  transformed  by 
(-l)n   before  transformation,  which  will  indeed  shift  the  spectrum  in 
such  a  way  that  its  origin  will  be  centered  in  the  transform  domain 
array.   In  this  investigation  a  different  approach  was  used  because 
Andrews'  method  still  leaves  an  unwanted  phase  shift  for  an  object  that 
is  centered  in  the  middle  of  the  array:  the  DFT  routine  asks  for  an 
object  that  is  centered  at  the  origin  for  an  output  without  phase  shift. 
The  center  of  the  array  was  defined  to  be  at  the  point  (l+N/2 ,  l+N/2). 
A  subroutine  (called  SHUFL,  see  Appendix)  was  written  which  interchanges 
the  four  quadrants  of  the  array  such  that  the  array  point  (l+N/2,  l+N/2) 
after  rearranging  occupies  the  point  (1,1).   (See  Fig.  4-7) 

In  interchanging  the  four  quadrants  of  the  array  by  a  simple 
translation  of  the  object,  advantage  was  taken  of  the  otherwise  undesira- 
ble postulated  existence  of  the  ghost  images.   After  transforming  the 
array,  the  transform  was  once  again  operated  on  by  SHUFL,  which  produced 
the  desired  effect  of  moving  the  origin  of  the  angular  spectrum  to  the 
center  of  the  array  (as  previously  defined)  without  any  unwanted 
associated  phase  shift. 
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B.   FAST  FOURIER  TRANSFORM 

The  discrete  Fourier  transform  requires  N2  complex  multiplications 
and  additions  for  a  one-dimensional  record  of  N  samples;  this  can  be 
prohibitive  for  even  moderately  large  N.   The  Fast  Fourier  Transform 
(FFT)  is  just  a  different  way  of  computing  a  DFT  (see  [Ref.  14],  but  as 
it  reduces  the  number  of  complex  multiplications  and  additions  on  the 
order  of  N*log2N  its  importance  can  hardly  be  overestimated.   The 
reduction  in  number  of  computations  is  even  more  dramatic  for  arrays 
of  higher  dimension.   Reference  [15]  lists  a  few  examples.   Without  the 
economy  of  the  Fast  Fourier  Transform  idea  the  DFT  would  be  a  useful 
tool,  but  to  be  used  very  sparingly  because  of  the  large  time  penalty 
to  be  paid.   Since  Cooley  and  Tukey  published  their  report  [Ref.  16] 
digital  signal  processing  has  received  new  impetus  because  transformation 
operations  are  now  possible  at  a  tolerable  even  though  still  not  a  very 
economical  time  expenditure.   The  FFT  algorithm  used  in  this  investiga- 
tion (subroutine  F0UR2)  [Ref.  17] ,  is  included  in  the  Appendix.   The 
time  needed  for  a  transformation  of  a  64  by  64  element  matrix  using 
F0UR2  was  less  than  2  seconds. 
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V.   THE  DETECTION  PROBLEM 

To  recapitulate,  the  problem  to  be  solved  was  to  reconstruct  the 
original  object  from  its  observed  sound  diffraction  pattern  without 
previous  knowledge  of  the  distance  of  the  object  from  the  observation  plane. 

A.   GENERAL  OBJECT  AND  UNIVERAL  RECONSTRUCTION  CRITERION 

For  an  unknown  object  distance  a  given  diffraction  pattern  can  be 
due  to  an  infinite  number  of 'objects'!   It  is  therefore  not  possible 
without  a  prior  knowledge  to  reconstruct  a  general  object  (not 
constrained  in  any  way  in  amplitude  and  phase)  by  using  a  simple 
criterion  function,  i.e.  one  which  does  not  contain  the  object 
distance  in  some  coded  form.   The  proof  of  this  statement  follows. 

Consider  an  object  U-^  belonging  to  a  rather  general  class  U 
of  objects,  all  of  which  are  related  to  each  other  by  some  common 
amplitude  and  phase  characteristics.   Suppose  also  there  is  a 
reconstruction  criterion  W  which  when  applied  to  all  diffraction 
patterns  A  from  all  objects  belonging  to  U  successfully  reconstructs 
the  respective  object  U^.   It  is  hoped  that  W  is  the  'universal 
reconstruction  criterion'. 

Suppose  a  diffraction  pattern  At  appears  at  a  distance  z  from  Ui  . 
Now  take  the  diffraction  pattern  A*  which  appears  at  a  distance  z*  in 
front  of  U,  to  be  a  new  object  to  be  reconstructed.   A*  gives  rise 
to  A.  at  a  distance  z,-z*  in  front  of  A*.   Application  of  W  on  A. 
will  always  reconstruct  U,  (according  to  the  second  hypothesis). 
Therefore  it  will  never  be  possible  to  extract  A*  from  A-^  by  using 
the  criterion  W.   So  W  cannot  be  the  desired  universal  reconstruction 
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criterion.   Since  U  and  W  were  chosen  arbitrarily  in  this  example, 
one  can  generalize  that  no  single  simple  criterion  can  be  used  to 
recover  all  objects. 

The  expression  'simple  criterion'  has  been  used  in  the  above 
discussion  for  the  reason  that  it  was  proposed  to  make  two  diffraction 
patterns  of  the  same  object  using  different  wavelengths.   Then,  on 
backward  propagation  both  images  should  coincide  at  the  object  location 
but  not  at  different  propagation  distances,  at  least  over  a  considerable 
portion  of  the  array.   This  idea  was  not  investigated  any  further 
because  it  was  thought  to  circumvent  the  issue  of  the  problem.   The 
difficulties  in  this  approach  are  obvious  at  once:  The  need  of 
recording  two  diffraction  patterns  with  the  associated  long  time 
expenditure,  the  need  of  almost  twice  as  much  computer  storage  and 
time  etc.   Moreover  it  was  hoped  that  a  single  diffraction  pattern  would 
contain  enough  information  to  find  the  original  object  if  the  class  of 
objects  to  be  reconstructed  were  somewhat  limited. 

Recently  a  reconstruction  time  has  been  proposed  by  Kalra  and 
Rodgers  [Ref.  18]  that  works  reliably  on  a  class  of  objects,  among  them 
the  point  scatterer  or  source,  which  have  a  magnitude  of  the  acoustical 
disturbance  which  reaches  its  maximum  at  the  object  location.   Unfortunately 
this  very  simple  criterion  will  not  reconstruct  successfully  objects 
whose  phase  distribution  departs  much  from  a  linear  one.   In  the  above 
reference  [Ref.  18]  the  possibility  of  extending  the  analysis  from  a  point 
source  to  a  general  object  is  suggested.   Of  course  any  object  can  be 
thought  of  as  a  (possibly  infinite)  sum  of  point  sources  and  the 
diffraction  pattern  due  to  that  object  as  resulting  from  the  appropriate 
accumulation  of  point  sources  if  the  phase  relationship  between  all 
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point  sources  is  accounted  for.   But  the  assumption  that  maximum 
amplitude  is  a  general  reconstruction  criterion  can  be  disproven 
(without  going  to  such  exotic  'objects'  as  other  diffraction  patterns) 
by  considering  an  uniformly  insonified  object,  whose  phase  varies  in 
a  quadratic  fashion  such  that  it  is  advanced  towards  the  edges  and 
retarded  in  the  center,  i.e.,  a  lens-type  phase  distribution.   Then  the 
sonic  energy  in  propagation  will  concentrate  near  the  focal  point  of 
this  'lens1.   As  the  energy  density  there  is  much  higher  than  at  the 
uniformly  insonified  object  (see  [Ref.  19]),  the  amplitude  of  the  sound 
field  must  be  higher  than  at  the  object.   This  is  due  to  the  constructive 
interference  introduced  by  the  chosen  phase  relationship.   The  situation 
as  above  was  simulated  on  the  computer:   A  16 A  by  8X   slit  with  a  phase 
distribution  that  had  an  equivalent  focal  length  of  10A  was  insonified 
by  a  unity  amplitude  plane  wave.   The  computed  amplitude  at  the  focal 
point  was  found  to  be  more  than  a  factor  of  8  higher  than  the  amplitude 
at  the  slit  plane.   So  the  proposed  'unbiased  estimator'  for  the  true 
object  distance  z*  of  above  reference  [Ref.  18]  will  always  choose  the  focal 
plane  as  object  location.   It  may  be  argued  that  a  lens-type  phase 
relationship  is  very  uncommon  in  nature,  but  random  phase  distributions 
(diffuse  scattering)  are  widely  observed,  and  random  phase  relationships 
between  the  point  sources  that  make  up  an  object  can  also  give  rise  to 
amplitudes  of  the  disturbance  which  are  higher  at  other  planes,  than 
at  the  object  location.   As  example  a  16A  by  8A  unity  amplitude  object 
was  taken  to  have  random  phase  distribution  between  array  points, 
uniformly  distributed  from  0  to  2tt  •   The  highest  amplitude  was  found 
to  be  2.1  times  the  value  at  the  slit  at  a  location  5A  in  front  of  the 
object.   This  must  not  necessarily  be  so;  for  other  distributions  the 
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amplitude  may  be  highest  at  the  object  location  (the  linear  distribution 
might  be  considered  a  very  special  case  of  a  random  distribution,  as 
might  be  the  parabolic  distribution) ,  but  it  introduces  an  element  of 
uncertainty  into  the  reconstruction  scheme  which  may  well  be  unbearable. 

In  the  initial  stages  of  this  investigation  it  was  hoped  that  cross- 
correlation  of  the  diffraction  pattern  with  a  point  source  at  the 
correct  distance  would  give  the  necessary  high  output  to  discriminate 
between  the  object  giving  rise  to  the  diffraction  pattern  and  some  other 
diffraction  pattern  due  to  this  object. 

Suppose  there  is  a  unity  amplitude  point  source  at  the  orgin.   Then 
Up  (x,y,o)  =  <S(x,y,o)   and 
Ap  (fx,fy,o)  =  1 
If  the  disturbance  due  to  this  point  source  is  propagated  a  distance  z 
then 

ApZ  (fx,fy,z)  =  Ap  (fx,  fy,  o)-  exp  [j^z  /l-Ufx)z-(Afy)z] 

=  1  •  exp  [j^J-z  /l-(Afx)MAfy)z] 

Now  cross-correlation  can  be  accomplished  by  multiplying  the  spectrum  of 

one  function  by  the  complex  conjugate  of  the  spectrum  of  the  other 

function  and  inverse  Fourier  transforming  the  result. 

Let  the  Fourier  transform  of  the  diffraction  pattern  U^Cx^jZ) 

be  given  by  A  (f  ,  f  ,  z) .   Then  the  result  of  the  proposed  crosscorre- 
^T.       X    y  re 

lation  can  be  written  as: 

R  (x.y.o)  =  F"1  {A  z  (fx,fy,z).  exp[-j^z  /l-(Afx)^-(Xfy) '-]} 
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which  is  recognized  as  the  result  obtained  by  backward  propagation  of 
the  diffracted  wavefront  by  a  distance  z,  namely  the  original  object. 
This  method  of  maximizing  the  crosscorrelation  cannot  be  used  as  it 
has  been  shown  that  the  amplitude  of  the  wavefront  at  the  object  need 
not  be  higher  than  everywhere  else  in  its  propagation  field. 

B.   OBJECTS  WITH  BOUNDARIES  AND  THEIR  RECONSTRUCTION 

In  the  search  for  a  class  of  objects  which  might  be  more  useful 

than  the  one  reconstructed  by  Kalra  and  Rodgers  [Ref.  18]  a  brief  study  of 

psychopictorics ,  that  subfield  of  psychophysics  which  deals  with 

natural  images  as  pictorial  stimuli  (see  Lipkin  and  Rosenfeld  [Ref.  20]) 

showed  that  the  three  most  important  psychopictorial  variables  are: 

contrast  and  border;  shape  and  geometry;  texture.   Of  these  contrast 

and  border  is  easiest  to  relate  to  what  is  actually  available  in  a 

digitized  picture.   To  treat  shape  and  geometry  would  mean  going  to  a 

pattern  recognition  system,  probably  by  using  a  crosscorrelation  scheme, 

with  the  attendant  necessity  of  storing  a  large  number  of  objects  to 

be  recognized.   Trying  to  base  a  detection  system  on  the  texture  of 

an  object  compared  to  that  of  its  surroundings  would  probably  be 

bound  to  fail  with  today's  still  rather  crude  acoustical  holography 

systems.   Thus,  trying  to  detect  objects  which  are  characterized  by 

contrast  and  boundary  seems  most  promising.   In  fact,  many  objects 

that  are  of  interest  to  acoustical  holography  fall  under  this  heading 

in  that  they  are  characterized  by  emitting  or  reflecting  large  amounts 

of  sonic  energy  and  being  separated  from  the  non-emitting  (or  less 

emitting)  background  by  a  rather  sharp  boundary  or  vice  versa. 

The  'picture'  of  a  mine  lieing  on  the  ocean  floor  will  show  a  sharp 

step  in  intensity  as  will  that  of  a  bone  imbedded  in  soft  tissue 

and  there  are  more  examples  that  might  be  cited. 
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1.  Objects  of  Interest 

It  is  hoped  then  that  this  class  of  objects  will  be  general 
enough  to  warrant  investigation.   Therefore  in  all  the  following 
discussion  it  will  be  assumed  that  the  object  plane  to  be  reconstructed 
from  its  diffraction  pattern  has  associated  with  it  (or  ideally  is 
characterized  by)  a  sharp  change  in  the  emitted  amplitude  (ideally  a 
discontinuity)  separating  regions  of  large  amplitude  from  those  of 
low  amplitude.   This  change  in  amplitude  will  not  be  assumed  to  be 
confined  to  a  single  point  or  to  a  number  of  disconnected  points 
because  the  word  'boundary'  or  'edge'  implies  connection  between  quite 
a  large  number  of  points.   In  this  sense  then  a  point  source  is  not 
a  valid  object  for  this  discussion.   A  line  object  cannot  be  hoped 
to  be  detected,  because  it  really  has  not  any  area  associated  with  it; 
the  same  objection  will  hold  for  an  object  that  is  checkerboard-shaped 
with  the  amplitude  changing  from  high  to  low  and  back  again  from  one  sampling 
point  to  the  next  one.   The  smallest  object  that  one  could  hope  to 
detect  is  composed  of  four  sampling  points  of  high  intensity  next 
to  each  other,  arranged  in  a  square.   This  probably  would  not  be 
a  great  obstacle  to  the  use  of  the  proposed  routine  as  it  implies  detection 
of  at  least  all  one  wavelength  by  one  wavelength  objects  if  the  Nyquist 
sampling  rate  is  used. 

2 .  Reconstruction  Schemes 

The  first  reconstruction  scheme  for  the  given  class  of  objects 
was  very  simple:  As  the  object  was  assumed  to  have  a  sharp  step  in 
amplitude  it  was  thought  that  the  difference  in  amplitude  between  one 
sampling  point  and  its  neighboring  sampling  point  would  have  a  maximum 
at  the  object  location  and  not  in  any  other  plane.   Therefore  that 
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plane  which  contains  this  maximum  would  be  taken  to  be  the  object 
plane  and  the  reconstructed  diffraction  pattern  in  this  plane  would  be  the 
image  of  the  original  object .   Unfortunately  the  above  assumption  was 
proven  to  be  not  valid  for  the  Nyquist  sampling  rate  if  the  phase 
distribution  over  the  object  departs  markedly  from  a  linear  one,  by 
assuming  such  an  object  and  simulating  its  diffraction  on  the  computer. 
Random  phase  objects  as  well  as  such  with  a  lens-type  phase  relationship 
showed  maximum  differences  between  some  sampling  points  which  were 
higher  (even  very  much  higher  in  the  case  of  the  parabolic  phase  object) 
in  other  planes  than  the  discontinuity  in  the  object  plane.   As  examples: 
A  16A  by  8X  unity-amplitude  object  with  uniformly  distributed  random 
phase  had  a  maximum  difference  of  1.73  at  a  distance  3  wavelengths  in  front 
of  the  object.   The  same  object  with  lens-type  phase  distribution  and 
a  focal  length  of  10X  showed  a  maximum  difference  of  3.85  in  the  focal 
plane.   At  first  this  was  not  taken  to  invalidate  the  whole  idea, 
because  theoretically  by  finer  and  finer  sampling  it  would  be  possible 
to  reduce  the  maximum  difference  between  sampling  points  in  any  other 
plane  so  far  that  it  would  be  smaller  than  the  given  discontinuity  in 
the  object  plane,  which  would  not  be  influenced  by  finer  sampling. 
But  apart  from  the  obvious  difficulty  in  an  actual  system  to  sample 
fractions  of  a  wavelength  apart  (the  detector  would  have  to  have  very 
small  area  indeed)  the  hoped  for  gain  could  only  be  realized  at  an  immense 
increase  in  amount  of  data  and  additional  computation  time  and  core 
memory  necessary  to  handle  that  many  data.   If  one  samples  finer  without 
increasing  the  number  of  samples  (which  might  be  possible  for  very 
small  objects)  nothing  is  gained.   This  hinges  on  the  fact  that  by 
sampling  above  the  Nyquiest  rate  too  few  frequency  components  will  be 
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available  to  faithfully  reconstruct  the  original  object.   (see  Sect.  IV) 
The  image  will  lose  its  clear  amplitude  step  and  instead  will  develop 
a  lower  amplitude  slope  and  overshoots;  in  this  way  the  characteristic 
that  discriminates  between  actual  object  and  its  diffraction  pattern  is 
lost.   Some  other  way  has  then  to  be  found  to  detect  the  edge  of  the 
intensity  distribution. 

In  the  above  scheme  no  use  was  made  of  the  fact  that  an  edge 
consists  of  more  than  one  spot  where  the  amplitude  changes  rapidly, 
that  a  number  of  such  changes  have  to  occur  at  adjacent  sampling 
points  and  that  these  places  of  rapid  change  have  to  separate  a 
region  of  high  amplitude  from  one  of  low  amplitude. 

The  problem  of  enhancing  or  detecting  an  edge  in  a  digital  picture 
has  been  treated  in  the  literature,  e.g. ,  by  Rosenfeld,  Lee,  and 
Thomas  in  [Ref.  21].   The  usual  treatment  is  to  take  differences  of  averages, 
e.g.  for  a  one-dimensional  case 


dmk  =  — 
m 


ak+m  +  •••+  ak+1  -  ak  -"•■  -ak-m+1 

where  a^  is  the  amplitude  of  the  point  being  investigated.   The 
difference  d,,=a,,,-a,  is  sensitive  to  every  local  fluctuation, 
whether  it  occurs  at  an  edge  or  not.   The  difference  d  .  for  large  m 
has  the  disadvantage  of  detecting  even  sharp  transitions  in  a  blurred 
manner  since  it  registers  a  difference  for  all  k  such  that  some  of  its 
terms  are  in  one  segment,  some  in  the  other  segment.   A  simple  method 
of  overcoming  these  difficulties  is  using  a  product  of  d  k's,  which 
Involves  both  small  and  large  m's.   This  product  can  only  be  large 
if  all  its  factors  are  large.   For  the  two-dimensional  case  this  edge 
detection  scheme  naturally  becomes  more  involved.   Operations  of  the  form 
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d  ■=    . !       J    (a        ,        +"*+a,       ,  +«««+a     , ,,       +*•• 

rs'nk  r(2s+l)        h+r'k+s  h+r  k-s  h+1  k+s 

+  ah+l*k-s)    ~    (ah'k+s+'"+ah'k-s4""+ah-r+l'k+s+"" 

+  V-r+l'k-s*' 
are  performed  and  products  of  the  results  for  various  large  and  small 

r's  and  s's  are  formed.   Again  the  product  can  only  be  relatively 
large  if  all  its  factors  are  large.   The  example  d   ,,,  shown  above 
is  used  for  the  detection  of  a  horizontal  edge,  the  parameter  s  increas- 
ing results  in  the  sensitivity  towards  orientations  departing  from  the 
horizontal  being  accentuated.   Each  possible  edge  orientation  has 
to  be  treated  by  a  different  combination  of  differences  of  averages, 
obviously  a  time-consuming  task  for  an  array  of  say  64  by  64  sample 
points.   The  edge  detection  operation  would  normally  then  be  followed 
by  a  curve  detection  routine  to  filter  out  any  spurious   output,  as  the 
output  of  the  edge  detection  step  will  be  noisy  for  an  edge  which  is  not 
perfectly  smooth  and  straight.   The  curve  detection  idea  is  based  on  the 
fact  that  the  values  resulting  from  the  edge  detection  should  be  high 
on  moving  along  the  edge  and  low  in  a  direction  perpendicular  to  it. 
In  other  words  a  high  point  on  an  edge  should  have  neighboring  high  points 
in  approximately  opposite  directions.   Points  which  do  not  satisfy  this 
criterion  are  regarded  as  spurious  and  discarded.   For  a  more  profound- 
discussion  of  this  topic  reference  is  made  to  Rosenfeld,  Thomas  and 
Lee  [Ref.  22].   Again,  the  time  needed  for  this  step  in  the  recognition 
process  cannot  be  regarded  as  negligible,  because  it  involves  checking  the 
neighborhood  of  N2  points  and  a  quite  involved  decision-making  process. 

To  receive  some  information  on  the  computation  time  needed  using 
the  above  method  of  edge  detection  a  first  approximation  to  it  was 
developed  and  tested.   In  order  to  make  the  operation  d   ,  ,  relatively 
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insensitive  to  edge  orientation  s  was  taken  to  be  1.   Differences 

d21'hk'd41'hk  and  d81'hk  were  formed  and  multiplied  together 
for  the  horizontal  and  vertical  edge  orientation  and  the  products 
for  the  horizontal  and  vertical  operation  added  together.   It  was 
hoped  in  this  way  to  overcome  the  necessity  of  investigating  all 
different  orientations.   The  result  obtained  with  this  technique  was 
promising.   (The  test  was  performed  using  a  circle,  whose  area  was 
filled  with  high  values,  as  a  boundary  comprising  all  edge  orientations 
as  input  to  the  routine.   Input  and  output  are  shown  in  Fig.  5-1  and 
Fig. 5-2).   Nevertheless  this  method  had  to  be  shelved  because  of  the  long 
computation  time  needed,  over  20  seconds  for  each  64  by  64  element 
array.   To  this  time  would  have  to  be  added  the  time  needed  for  the 
previously  discussed  curve  detection  algorithm  also. 

The  difficulty  of  having  to  treat  each  orientation  separately 
is  analogous  to  that  which  is  encountered  in  trying  to  accomplish  edge 
detection  by  a  cross-correlation  scheme,  which  also  is  very  sensitive 
to  a  relative  rotation  between  the  reference  object  and  the  object  to 
be  tested,  and  additionally  has  to  account  for  different  phase  relation- 
ships also. 

In  order  to  surmount  these  difficulties  due  to  various  edge  orienta- 
tions a  new,  rather  simple  (and  thus  more  error-prone)  procedure  was 
worked  out.   In  this  algorithm  first  the  largest  difference  between  any 
two  adjacent  array  points  in  the  plane  under  investigation  is  found. 
This  is  used  to  set  two  thresholds,  an  upper  (THRUP)  and  a  lower  one 
(THRLO) .   These  thresholds  are  kept  variable,  so  that  they  can  be  adjusted 
according  to  the  kind  of  object  one  is  looking  for  (if  such  previous 
information  is  available).   A  setting  of  60%  of  the  maximum  difference 
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Figure  5-1 
Amplitude  Contour  Plot  of  Input  to 
Edge  Detection  Scheme  to  be  Tested 


64 


Figure  5-2 
Amplitude  Contour  Plot  of  Output  of 
Edge  Detection  Scheme  to  be  Tested 
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between  array  apoints  for  the  upper  threshold  and  20%  for  the  lower 
one  seems  in  general  to  give  good  results.   An  edge  in  one  direction 
is  assumed  to  exist  at  an;  array  point  (i,k)  if  for  a  downgoing  step 

(a^j^-a. ,.,-,)  is  greater  than  the  upper  threshold  and  both  | a. .  av+l- a± »V+2  1 

and  la     -a    I  are  less  than  the  lower  threshold.   This  criterion  is 
1  i ' k-1  i ' k ' 

analogous  to  the  edge  detection  routine  described  earlier  with  S  =  0  and 
the  product  operation  replaced  by  its  logical  equivalent  'AND' :  The 
output  in  both  cases  will  be  positive  ('YES'  in  this  case,  large  in  the 
former)  if  all  the  inputs  (differences)  to  the  operation  are  positive 
('YES'  or  large) . 

This  edge  criterion  is  used  twice  at  each  array  point,  for  horizontal 
and  for  vertical  differences.   In  this  manner  corners  as  points  belonging 
to  two  edges  at  the  same  time  are  weighted  more  than  simple  edge  points. 
(The  deemphasis  of  corners  is  one  unfortunate  feature  of  the  earlier  edge 
detection  scheme).   Also,  a  distinction  is  made  between  upward  steps 
((a.  ,,  , ,  -a.,,)  must  be  greater  than  the  upper  threshold,  while  the  other 
two  conditions  stay  unchanged  in  this  case)  and  downward  steps.   This 
feature  is  used  in  that  part  of  this  scheme  which  assigns  different  values 
to  each  edge  point  according  to  whether  it  is  connected  to  other  edge 
points  or  not:   The  step  heights  for  an  array  point  (i,k)  are  computed  in 
both  the  horizontal  and  vertical  directions;  values  of  positive  going 
steps  are  kept  separate  from  negative  steps;  then  the  four  array  points 
next  to  (i,k)  which  have  been  treated  previously  in  the  sequential  edge 
search  are  tested  for  the  highest  absolute  value  of  assigned  'edgemeasure' 
If  this  is  due  to  an  upward  going  edge  the  height  of  the  positive  step  at 
the  point  in  question  is  added  to  it  and  stored  at  (i,k).   The  treatment 
is  analogous  if  the  highest  value  stored  next  to  (i,k)  has  been  found 
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from  a  negative  edge;  the  negative  step  height  at  (i,k)  is  added  to  it 
and  kept  at  (i,k).   If  the  highest  adjoining  edge  measure  is  positive, 
but  only  negative  steps  occur  at  (i,k)  this  negative  step  height  is 
stored  at  (i,k)  without  any  addition  of  previously  computed  edge  measure. 
The  same  thing  happens  if  the  adjoining  array  points  contain  zero. 
Zero  is  stored  at  (i,k)  regardless  of  adjoining  edge  measures  if  the 
situation  at  (i,k)  and  its  vicinity  does  not  meet  the  edge  criterion 
using  the  upper  and  lower  thresholds.   For  a  more  exhaustive  treatment 
of  all  different  cases  that  may  arise  see  the  computer  flow  diagram  for 
subroutine  ADERIV  in  the  Appendix.   After  the  whole  array  has  been 
treated  in  this  manner,  all  edge  measures  computed  are  added  in  absolute 
value;  the  result  is  taken  to  be  the  overall  edge  measure  of  the  given 
array. 

It  must  be  admitted  that  the  only  justification  for  this  scheme 
that  can  be  given,  is  that  it  seems  natural  in  looking  for  edges  to 
discriminate  between  large  and  small  differences  between  adjoining 
points  and  that  longer  connected  boundaries  should  be  weighted  more  than 
short  ones:   three  points  of  rapid  change  next  to  each  other  might  be 
considered  at  least  as  a  beginning  of  an  edge,  whereas  two  such  points 
could  be  just  a  random  occurrence.   The  normal  height  of  the  thresholds 
as  given  is  found  entirely  from  practice  with  this  routine  and  has  no 
profound  mathematical  background.   It  may  be  argued  that  in  the 
decribed  algorithm  too  much  importance  is  attributed  to  connection 
between  edge  points  and  too  little  to  the  height  of  the  single,  steps . 
But  it  has  been  observed  that  if  an  array  containing  fewer,  large  steps 
is  compared  with  another  array  with  many  smaller  steps  next  to  each 
other,  the  latter  is  more  likely  to  b e  the  looked  for  object  than  the 
former.   Besides  that,  a  shift  in  the  relative  weight  balance  between 
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'connection'  and  'step  height'  is  easily  accomplished  by  multiplying 
the  largest  absolute  value  of  edge  measure  adjoining  the  point  (i,k)  by 
some  factor  less  than  one  before  the  addition  and  storage  at(i,k). 

The  algorithm  can  easily  be  made  more  discriminating  against  steps 
which  do  not  belong  to  an  actual  edge  by  taking  into  account  not  merely 
two  differences  next  to  the  point  in  question  but  investigating  more 
points  in  its  neighborhood.   Then  it  must  be  borne  in  mind  that  as  more 
small  differences  in  one  direction  are  required  to  satisfy  the  step 
criterion  the  smallest  detectable  object  will  become  larger.   Another 
means  of  increasing  the  rejection  capability  is  by  increasing  the 
upper  threshold  and  reducing  the  lower  threshold  if  that  is  compatible 
with  the  kind  of  object  to  be  investigated:  setting  a  lower  threshold 
of  0.1  of  the  maximum  difference,  if  the  amplitude  of  the  background  of 
the  object  is  uniformly  random  distributed  between  zero  and  20%  of  the 
object  discontinuity,  could  reject  quite  a  number  of  bona  fide  edge  points 
thus  reducing  the  probability  of  detection.   It  is  inadvisable  to  work 
with  an  edge  detection  routine  which  is  too  discriminating  from  another 
point  of  view  also:  this  routine  might  not  pass  the  diffraction  patterns 
a  small  distance  behind'  and  in  front  of  the  object  plane  with  a  high 
enough  edge  measure  to  make  them  obvious.   Thus  in  stepping  through  the 
volume  of  interest  the  computer  would  have  to  investigate  very  near  to 
the  actual  object  to  receive  an  indication  of  its  presence  or  it  would 
be  missed  altogether  and  the  computer  would  settle  for  some  other  than 
the  true  object  distance  and  output  the  diffraction  pattern  there  as 
the  reconstructed  object.   This  would  mean  that  many  small  steps  would 
have  to  be  taken,  which  implies  very  long  computation  times.   This  case 
could  not  appear  if  the  overall  edge  measure  as  computed  above  were 
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steadily  increasing  towards  the  actual  object  location.   Then  it  would 
be  just  a  matter  of  finding  the  maximum  edge  measure  and  the  object 
would  be  detected.   Unfortunately,  the  situation  is  not  that  simple. 
Figure  5-3  shows  the  computed  edge  measures  for  a  16  by  8  wavelength, 
unity  amplitude  object  with  random  phase,  versus  propagation  distance 
in  steps  of  2.5  wavelengths.   It  shows  quite  a  number  of  planes  in 
which  the  edge  measure  has  a  relative  maximum.   This  is  not  only  due 
to  the  way  in  which  the  edge  measure  is  computed  but  also  due  to  the 
rapid  variation  of  the  diffracted  wavefront  with  distance. 

A  reconstruction  method  then  has  to  s  tart  from  a  diffraction  pattern 
and  taking  advantage  of  all  previous  information  (which  might  have  been 
obtained  by  rangegating)  propagate  this  diffracted  wavefront  back  to  some 
starting  distance.   From  there  to  some  given  final  distance  the  wavefront 
will  be  propagated  back  in  as  small  a  number  o f  steps  (as  large  a  step 
size)  as  possible  without  missing  the  object  altogether,  meaning  that 
the  step  size  has  to  be  so  small  that  an  indication  of  a  possible  object 
is  received  by  a  relative  edgemeasure  maximum  less  than  half  a  step  size 
behind  or  in  front  of  the  actual  object  location.   In  the  proposed 
algorithm  this  is  accomplished  by  setting  the  parameters  ZFIRST  for  the 
starting  distance,  ZLAST  for  the  final  distance  and  NRSTEP  for  the 
number  of  steps  to  be  used.   The  choice  of  step  size  for  this  first  pass 
through  the  volume  of  interest  is  thus  left  to  the  operator.   For  an 
object  of  16  by  8  wavelengths  and  Nyquist  rate  sampling  the  step  size 
should  not  be  chosen  larger  than  about  2  wavelengths.   Larger  objects 
can  bear  a  larger  step  size. 
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Figure  5-3 
Edge  Measure  of  Diffraction  Pattern  of  Random 
Phase  Slit  Versus  Propagation  Distance 
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All  relative  maxima  in  overall  edge  measure  observed  on  this  first 
pass  have  to  be  stored  for  future  reference,  together  with  the  propaga- 
tion distances  at  which  they  occur. 

Different  treatments  of  the  result  of  this  first  pass  are  possible 
now.  For  instance  the  reconstructed  wavefront  at  all  locations  of  a 
relative  maximum  might  be  displayed  on  a  plotter,  or  graphics  terminal 
and  the  decision  as  to  the  actual  location(s)  of  object(s)  left  to  the 
experience  of  the  operator,  who  would  then  investigate  those  locations 
closer.  The  method  chosen  for  the  proposed  algorithm  differs  from  the 
above  by  making  further  use  of  the  computer's  decision-making 
capabilities. 

After  the  first  pass  through  the  whole  distance  of  interest  the 
vector  of  edge  measure  maxima  is  reordered  according  to  their  size  and 
those  possible  object  locations  discarded,  whose  edge  measure  is  only  10% 
or  less  of  the  maximum  edge  measure  encountered  on  the  first  pass,  it 
being  assumed  that  closer  investigation  of  these  locations  will  not 
yield  the  looked  for  object.   The  remaining  locations  are  investigated 
in  a  fixed  step  size  which  is  set  equal  to  the  sampling  increment,  i.e, 
to  one  half  a  wavelength  if  the  Nyquist  rate  is  used,  from  a  distance 
one  half  of  the  original  step  size  in  front  of  the  relative  maximum  to 
one  half  that  step  size  behind  this  location.   The  absolute  maximum  in 
edge  measure  in  this  interval  is  found  and  stored  in  a  different  array. 
After  all  relative  maxima  resulting  from  the  first  pass  have  been 
investigated  the  array  of  second  pass  maxima  is  scanned  for  the  overall 
maximum  and  those  locations  discarded  where  the  edge  value  is  below 
a  certain  fraction  of  the  overall  maximum.   This  is  accomplished  by 
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setting  a  threshold  THRCON  to  a  fraction  of  1.0:  one  would  not  expect 
maximum  of  less  than  0.2  to  0.3  of  the  overall  maximum  to  be  due  to 
actual  objects.   Also  incorporated  into  the  algorithm  was  the  possibi- 
lity of  investigating  fewer  relative  maxima  resulting  from  the  first 
pass  by  reading  in  a  value  for  NRMAX,  e.g.,  if  it  is  known  that  the 
object  exists  only  in  a  single  plane  parallel  to  the  observation  plane. 
This  option  should  be  used  with  caution  because  there  is  no  guarantee 
that  the  highest  relative  maximum  after  the  first  pass  will  be  associated 
with  the  actual  object  for  all  objects  (although  this  is  true  for  many 
simple  objects). 

The  results  of  propagating  the  wavefront  to  the  locations  of  the 
remaining  second  pass  maxima  are  then  plotted  on  the  CALCOMP  -  plotter, 
either  in  form  of  a  contour  plot  of  the  computed  amplitude  or  as  a 
surface  composed  of  the  values  of  amplitude  at  all  array  points,  from 
which  the  hidden  lines  are  removed  pother  graphical  output  devices  can 
also  be  used).   The  contour  plot  will  give  a  much  better  pictorial 
representation  of  the  outline  of  the  object,  but  a  good  part  of  the 
amplitude  information  is  lost  by  plotting  only  a  number  of  contour  lines. 
The  three-dimensional  representation  is  t  o  be  preferred  on  the  basis 
of  information  contained  in  the  plot  (it  also  suffers  from  loss  of 
information  because  it  is  not  possible  to  see  through  the  surface) 
but  it  takes  much  longer  computation  times  and  requires  much  more 
core  memory  than  the  contour  plot.   Also  it  must  be  remembered  that 
the  apparent  dimensionality  is  not  due  to  any  three-dimensional  object 
of  that  shape  but  rather  represents  the  relative  amplitude  of  the 
wavefront  in  a  single  plane.   Subroutines  used  for  this  output  are  for 
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the  contour  plot  CONTUR  [Ref.  23]  and  PR01  for  the  surface  plot 
[Ref.  24].   For  a  64  by  64  element  array  to  be  treated  CONTUR  takes 
on  the  order  of  50  sec  if  six  contour  lines  are  plotted;  PR01  needs 
for  the  same  array  from  50  sec  up  to  10  min  depending  on  how  complicated 
the  surface  to  be  plotted  is. 
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VI.   TESTS  AND  RESULTS 

Part  of  the  tests  performed  were  designed  to  check  the  propagation 
algorithm,  the  greater  part  were  used  to  try  the  automatic  reconstruction 
scheme.   For  this  latter  purpose  cumputer  generated  diffraction  patterns 
and  one  set  of  actually  observed  data  were  used. 

A.    TEST  OF  THE  PROPAGATION  ALGORITHM 

As  the  spatial  frequency  approach  to  diffraction  is  limited  in  the 
maximum  useful  propagation  distance  due  to  the  presence  of  the  postulated 
ghost  images  it  is  not  possible  to  apply  the  usual  test  for  diffraction 
algorithms  in  which  the  amplitude  distribution  at  a  very  large  progaga- 
tion  distance  has  the  functional  form  of  a  Fourier  transform  of  the 
diffracting  aperture.   Instead  the  fact  was  used  that  illuminating  an 
aperture  with  a  spherical  wave  converging  to  a  point  PQ  will  give  rise 
to  an  amplitude  distribution  in  a  plane  perpendicular  to  the  direction 
of  propagation  and  containing  P0,  which  again  is  simply  related  to 
the  Fourier  transform  of  the  diffracting  aperture.   [Ref.  7]. 

A  10X  by  10X  square  aperture  was  simulated  to  be  insonified  by 
a  spherical  wave  converging  to  a  point  40A  behind  the  aperture.   Using 
the  proposed  propagation  algorithm  and  propagating  the  diffracted 
wavefront  40A  resulted  in  the  output  of  Fig.  6-1.   It  is  readily 
apparent  that  this  output  has  the  functional  form  of  (sin(x)/x) • (sin(y)/y) 
(which  is  the  Fourier  transform  of  a  uniformly  insonified  rectangular 
aperture).   The  fact  that  the  surface  of  Fig.  6-1  is  not  very  smooth  is 
due  to  the  relatively  small  number  of  surface  points.   Due  to  the 
presence  of  the  ghost  images  the  output  is  deteriorated  towards  the 
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boundary  of  the  presented  surface  (the  amplitude  does  not  fall  off 
as  fast  as  it  should  in  following  a  (sin(x)/x) • (sin(y)/y)  functional 
relationship) . 

It  is  of  interest  to  note  that  even  after  about  30  sequential 
propagation  steps  the  result  differs  from  that  obtained  in  one  step 
only  in  the  fourth  decimal  place.   This  good  agreement  could  never  be 
expected  if  the  Fresnel  equation  were  used. 

B.    TESTS  OF  THE  RECONSTRUCTION  SCHEME 

1.   Tests  Using  Computer  Generated  Diffraction  Patterns 

Most  of  the  tests  of  the  automatic  reconstruction  scheme  were 
performed  using  computer  generated  diffraction  patterns.   They  have 
the  advantage  of  being  perfect  in  the  sense  of  having  no  distortion 
introduced  by  noise  in  the  detection  system,  by  reverberation  and  by 
scan  non-uniformities.   In  all  tests  involving  computer  generated 
diffraction  patterns  the  distance  from  object  to  observation  plane 
was  chosen  arbitrarily  (subject  to  the  constraint  of  being  small 
enough  for  the  results  of  the  spatial  frequency  approach  to  diffraction 
to  be  valid).   The  numbers  of  steps  given  in  the  following  discussions 
are  those  for  the  first  pass  through  the  volume  of  interest  (i.e, 
those  given  by  setting  the  input  NRSTEP) .   The  number  of  steps  for  the 
closer  investigation  of  relative  maxima  in  edge  measure  is  computed  by 
the  detection  algorithm  (as  discussed  in  the  last  section) .   All  com- 
putation times  are  without  those  necessary  for  the  output  subroutines 
(CONTUR  or  PR01)  . 

Case  a:   The  simplest  object  whose  reconstruction  was  attempted  was  a 
16X  by  8X  uniformly  insonified  rectangular  slit  with  all  its 
sample  points  in  phase  thus  simulating  a  specular  object.   The  diffraction 
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pattern  was  observed  175. 3X  in  front  of  the  object;  in  the  recon- 
struction the  closest  distance  to  be  investigated  was  set  to  160X, 
the  furthest  distance  was  192X;  this  volume  of  interest  was  investi- 
gated in  12  steps.   A  reconstruction  (see  Fig.  6-2)  was  obtained  at 
a  distance  of  175.25  wavelengths  from  the  observation  plane, 
i.e.  only  0.05X  away  form  the  original  object.   The  time  needed  was 
about  2.6  min. 

Case  b:  The  same  object  as  in  case  a  was  assumed  with  the  additional 
difficulty  of  a  background  insonification  with  an  amplitude  which  was 
randomly  distributed  between  zero  and  0.2  of  the  object  amplitude. 
Also  the  phase  of  the  background  insonf ication  was  taken  to  be 
uniformly  random  distributed  between  zero  and  2tt,   The  distance 
between  object  and  observed  diffraction  pattern  was  177. 2X.   The 
volume  of  interest  and  number  of  steps  were  as  in  case  a.   The  obtained 
reconstruction  of  Fig.  6-3  was  located  0.05X  behind  the  actual 
object  location.   In  order  to  show  the  random  amplitude  variation  of  the 
background  here  the  output  is  given  in  form  of  a  3-dimensional  surface. 
Computation  time  was  about  2.6  min. 

Case  c:   An  object  of  the  same  size  and  shape  as  in  case  b,  was  treated 
in  this  case.   The  random  background,  too,  was  the  same  as  in  case  b. 
This  case  differed  from  case  b,  in  the  assumption  of  a  parabolic  phase 
relationships  for  the  object  which  gave  it  an  equivalent  focal  length 
of  10  wavelengths.   The  observation  plane  was  located  177. IX  in  front 
of  the  object.   Initial  and  final  distances  remained  unchanged  from  cases 
a  and  b,  as  did  the  number  of  steps.   An  object  reconstruction  was 
obtained  at  a  distance  of  177.25,  requiring  a  computation  time  of  3.3 
min. 
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Figure  6-2 
Amplitude  Contour  Plot  of 
Reconstructed  Specular  Object 
(Case  a) 
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Case  d:   To  show  that  the  proposed  reconstruction  scheme  works  also  for 
other  than  rectangular  objects  a  shape  was  chosen  for  the  input  object 
which  can  best  be  visualized  by  viewing  the  reconstructed  output  of 
Fig.  6-4.   Background  and  phase  distribution  of  the  object  are  the  same 
as  in  case  c.   Propagation  distance  to  the  observation  plane  was  177. IX. 
Initial  and  final  distances  for  the  investigation  were  set  to  160.0  and 
190.  .0  wavelengths  respectively,  the  investigation  was  performed  in  12 
steps.   The  reconstructed  object  (see  Fig.  6-4)  was  located  177. 25A 
from  the  observation  plane.   Time  needed  for  the  reconstruction  was 
about  4  min.   Figure  6-5  shows  the  same  output  as  Fig.  6-4,  only  in  the 
3-dimensional  surface  plot.   This  figure  was  included  here  to  illustrate 
some  of  the  shortcomings  of  this  output  scheme:   Even  for  a  relatively 
simple  shape  as  that  from  Fig.  6-4  the  3-dimensional  plot  does  not 
show  the  object  shape  clearly.   Moreover  it  exhibits  certain  errors 
due  to  the  rapid  variation  of  the  amplitude  function  which  cannot 
sufficiently  be  accounted  for  in  the  interpolation  scheme.   The  contour 
plot  (Fig.  6-4)  on  the  other  hand,  does  not  show  the  random  variation 
of  amplitude  in  the  background  insonif ication  while  exhibiting  quite 
clearly  the  shape  of  the  object. 

Case  e:   For  this  case  a  16A  by  8A  rectangular  uniformly  emitting 
object  was  assumed  to  have  a  uniformly  distributed  random  phase  distri- 
bution between  zero  and  2tt  whcih  simulates  a  diffuse  object.   No  back- 
ground sound  was  taken  to  exist.   The  diffraction  pattern  was  calculated 
for  a  observation  distance  of  177. 5A.   The  volume  of  interest  remained 
unchanged  between  160A  and  190A,  as  did  the  number  of  steps  in  the 
investigation.   One  would  expect  this  case  to  be  harder  to  solve  due 
to  the  strong  deterioration  of  the  best  possible  reconstructed  object 
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Figure  6-4 
Amplitude  Contour  Plot  of  Reconstructed  Parabolic  Phase 
Object  with  Random  Phase,  Random  Amplitude  Background  (Case  d) 
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(see  Fig.  3-6).   The  computer  found  two  possible  object  locations. 
Plots  for  the  diffraction  patterns  at  these  locations  (177. 5\   and 
188. 25A)  are  shown  in  Fig.  6-6  and  6-7  respectively.   The  first  of  these 
locations  is  the  actual  object  location  and  the  plot  Fig.  6-6  the 
best  output  that  can  be  expected.   Computation  time  was  5.5  min.   It 
is  surprising  that  the  edge-detection  routine  does  not  discard  the 
array  of  Fig.  6-7,  but  instead  assigns  it  a  value  of  overall  edge- 
measure  which  is  70%  of  that  of  the  actual  object.   Here  it  would  be 
necessary  for  a  human  operator  to  use  his  pattern-recognition  capability 
to  discard  the  second  object  hypothesis. 

Case  f :  As  the  last  test  for  the  reconstruction  of  objects  which  exist 
only  in  a  single  plane  a  16A  by  8X   object  was  simulated  with  a 
background  insonif ication  whose  amplitude  was  uniformly  random  distributed 
between  zero  and  0.2  and  whose  phase  was  assumed  to  vary  uniformly  between 
zero  and  2tt.   The  object  itself  had  an  amplitude  which  consisted  of 
a  constant  term  of  0.9  to  which  was  added  a  randomly  varying  value 
which  was  uniformly  distributed  between  zero  and  0.2.  Also,  a  random 
phase  variation  between  zero  and  2tt  was  assumed  for  the  object.   The  actual 
distance  between  object  and  observation  plane  was  177. 5X.   The  initial 
distance,  final  distance  and  number  of  steps  were  left  the  same  as  in 
case  c.   The  computer  required  5  min.  to  calculate  two  tentative  object 
locations  at  177. 5A  and  163. 75A.   The  diffraction  patterns  at  these 
locations  are  shown  in  Figures  6-8  and  6-9  respectively.   The  edge- 
measure  assigned  by  the  computer  to  the  spurious  output  of  Fig.  6-9  was 
0.75  of  that  of  the  correct  output. 

From  all  these  tests  it  would  seem  that  one  could  expect  the 
computer  to  be  able  to  reconstruct  any  object  which  has  a  sharp  edge 
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Figure  6-6 
Amplitude  Contour  Plot  of  Reconstructed  Random  Phase  Object 

(Case  e) 


84 


Figure  6-7 
Amplitude  Contour  Plot  of  Spurious  Output 
(Case  e) 
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Figure  6-8 
Amplitude  Contour  Plot  of  Reconstructed  Random  Phase  Object 
with  Random  Amplitude  Variation.   Background  Assumed  to  have 
Random  Phase,  Random  Amplitude  (Case  f) 
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Figure  6-9 
Amplitude  Contour  Plot  of  Spurious  Output 

(Case  f) 
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and  exists  only  in  a  single  plane  with  rather  high  probability. 
Naturally  each  output  would  have  to  be  scrutinized  by  an  operator 
such  that  spurious  output  which  could  appear  in  addition  to  the  actual 
object  reconstruction  would  be  discarded.   The  proposed  edge  detection 
routine  is  not  infallible,  but  it  performs  its  operation  in  rather 
short  time  which  was  the  main  reason  for  setting  it  up  in  the  existing 
manner . 

The  next  set  of  tests  performed  consisted  of  trying  to  reconstruct 
input  objects  existing  in  two  different  parallel  planes.   As  in  backward 
propagation  the  real  image  is  reconstructed,  the  first  difficulty  is 
immediately  obvious:  the  diffraction  pattern  converges  to  that  existent 
in  the  plane  of  the  object  on  backward  propagation,  but  the  sound  waves 
emanating  from  the  object  do  not  stop  on  the  reconstructed  object  but 
are  propagated  farther  back  on  the  next  propagation  step.   This  poses 
no  great  problem  for  an  object  that  exists  only  in  a  single  plane  parallel 
to  the  observation  plane.   But  consider  an  object  consisting  of  two 
extended  sound  sources  in  two  different  parallel  planes.   The  sound  from 
the  source  farther  from  the  observation  plane  is  propagated  into  the 
plane  of  the  nearer  source  there  to  be  absorbed,  phase  shifted  and/or 
diffracted,  as  the  case  may  be.   The  sound  field  that  then  exists  in  the 
plane  of  the  front  source  is  propagated  to  the  observation  plane.   On 
backward  propagation  the  field  converges  to  that  previously  existent 
in  the  front  plane.   So  far,  no  fault  can  be  found  with  the  described 
situation.   But  on  further  backward  propagation  (the  second  source  should 
also  be  reconstructed)  the  sound  field  introduced  by  the  front  source 
is  also  propagated  back  into  the  plane  of  the  farther  object,  there  to 
give  rise  to  a  disturbance  which  did  not  exist  in  actuality.   In  fact 
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this  is  not  too  different  from  the  case  of  the  optical  hologram  and  the 
eye  because  the  eye  converts  the  virtual  image  to  a  real  image  on  the 
retina.   By  focusing  on  a  part  of  the  image  that  is  in  a  plane  farther 
away,  that  part  which  is  nearer  is  still  seen,  in  a  blurred  fashion, 
in  that  same  plane,  but  the  brain  rules  out  the  hypothesis  of  a 
composite  sharp  and  blurred  object:   The  unwanted  information  is  largely 
ignored.   How  recognition  is  accomplished,  i.e.,  why  certain  object 
hypotheses  are  accepted,  other  not,  remains  unknown  to  date,  even 
though  the  necessity  of  experience  is  generally  assumed  (apart  from 
the  obvious  advantage  of  stereoscopic  vision  if  the  distance  is  not  too 
large).   Huang  [Ref.  25]  and  Ichioka,  Izumi  and  Suzuki  [Ref.  26]  show  a 
method  of  generating  diffraction  patterns  due  to  such  objects.   The 
disturbance  due  to  the  back  object  is  propagated  into  the  plane  of  the 
front  object.   There  that  part  of  the  resulting  disturbance  that 
corresponds  to  the  extension  of  the  front  object  is  replaced  by  the 
disturbance  due  to  the  front  object  and  the  resulting  composite 
disturbance  is  propagated  into  the  observation  region.   This  scheme 
was  used  for  generating  the  necessary  diffraction  patterns  in  this 
investigation. 

Case  g:   The  simplest  object  of  this  category  consisted  of  two  16X 
by  8X  rectangles,  assumed  to  be  specularly  emitting  (i.e.,  all  sample 
points  were  in  phase).   These  objects  were  located  in  planes  which  were 
15\  behind  one  another.   One  rectangle  was  rotated  90°  with  respect  to 
the  other  around  the  propagation  axis  and  displaced  in  such  a  manner 
that  the  projections  of  both  rectangles  along  the  axis  of  propagation 
overlapped  by  only  a  small  amount  (about  0.1  of  the  area  of  one 
rectangle).   The  actual  object  configuration  can  very  well  be 
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visualized  by  studying  Figs.  6-10  and  6-11,  which  are  the  results  of 
the  computer  reconstruction.   The  diffraction  pattern  was  observed 
100.4  wavelengths  in  front  of  the  front  rectangle.   The  volume  of 
interest  was  taken  to  extend  from  85X  to  130X  from  the  diffraction 
pattern.   Therefore  the  number  of  steps  in  the  investigation  had  to  be 
increased  to  30.   The  reconstructions  obtained  were  located  at  100. 25X 
and  115. 25X  from  the  observation  plane,  i.e.,  0.15X  in  front  of  the 
first  and  second  rectangle  respectively.   The  total  computation  time 
was  about  8.5  min. 

Case  h:   The  same  geometry  as  for  case  g  was  assumed  but  here  both 
rectangles  had  a  random  phase  distribution  uniformly  distributed  from 
zero  to  2ir.   The  volume  of  interest  and  number  of  steps  for  the 
investigation  were  not  changed  from  case  g,  but  the  observation  plane 
was  located  99.63X  in  front  of  the  front  rectangle.   The  computer 
obtained  8  object  locations  (at  93.53X,  99.69X,  103. 85X,  106. 95X, 
116. 26X,  120. 91X,  123. 79X  and  126. 12X),  a  clear  indication  that  this 
reconstruction  was  not  quite  successful.   The  array  corresponding  to 
the  actual  location  of  the  front  rectangle,  Fig.  6-12,  was  only  third  in 
overall  edge-measure.   An  output  which  is  typical  for  the  spurious 
objects  computed  is  shown  in  Fig.  6-13.   The  amplitude  distribution  of 
Fig.  6-13  was  the  fourth  highest  in  edgemeasure;  it  is  computed  at  a 
distance  of  7.32X  behind  the  front  object.   Figure  6-14  shows  the  output 
which  is  the  lowest  of  all  eight  computed  in  edgemeasure.   Its  generating 
diffraction  pattern  is  located  1.6X  behind  the  second  rectangle.   It  is 
not  evident  that  this  output  is  strongly  related  to  the  second  rectangle 
(a  comparison  with  Fig.  6-11  is  helpful).   The  existence  of  the  front 
rectangle  can  be  seen  quite  clearly  in  Fig.  6-12  (compare  with  Fig.  6-10). 
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Figure  6-10 
Amplitude  Contour  Plot  of  Reconstructed  Diffraction 
Pattern  in  Front  Plane  of  Double  Specular  Object 

(Case  g) 


91 


Figure  6-11 
Amplitude  Contour  Plot  of  Reconstructed  Diffraction 
Pattern  in  Back  Plane  of  Double  Specular  Object 
(Case  g) 
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Figure  6-12 
Amplitude  Contour  Plot  of  Reconstructed  Diffraction 
Pattern  in  Front  Plane  of  Double  Random  Phase  Object 

(Case  h) 
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Figure  6-13 
Amplitude  Contour  Plot  of  Spurious  Output 

(Case  h) 
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Figure  6-14 
Amplitude  Contour  Plot  of  Reconstructed  Diffraction  Pattern 
Near  Back  Plane  of  Double  Random  Phase  Object 

(Case  h) 
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On  studying  Fig.  6-12  and  6-13  it  is  not  surprising  that  the  computer 
cannot  clearly  recognize  the  output  of  Fig.  6-12  as  the  actual  object. 
The  deterioration  of  the  output  due  to  the  loss  of  evanescent  waves  is 
very  apparent  in  a  situation  as  complicated  as  the  one  treated  in  this 
case. 

Case  i:   In  order  to  investigate  the  effects  of  varying  phases  between 
both  rectangular  objects  and  those  of  larger  amounts  of  overlap,  two 
rectangles  of  unity  emission  amplitude  were  simulated.   These  were 
assumed  to  be  emitting  in  phase,  but  their  separation  in  the  direction 
of  propagation  was  taken  to  be  an  odd  number  of  half  wavelengths, 
namely  14. 5X.   Both  rectangles  were  of  size  16X  by  8X,  one  was  rotated 
90°  versus  the  other  around  the  propagation  direction  and  both  were 
centered  in  the  input  array.   The  overlap  thus  became  one  half  the  area 
of  one  rectangle.   The  diffraction  pattern  was  observed  at  100. AX  from 
the  front  object.   Investigation  took  place  in  30  steps  between  85X  and 
130X  from  the  observation  plane.   The  computer  generated  two  possible 
object  locations  at  100. 75X  (Fig.  6-15)  and  at  113. 5X  (Fig.  6-16)  from 
the  diffraction  pattern.   This  took  6.2  min.  of  computation  time.   The 
front  rectangle  being  reconstructed  0.35X  behind  its  actual  location 
shows  up  very  well.   Due  to  the  destructive  interference  between  front 
and  back  object  the  location  of  the  second  rectangle  could  not  be  found 
exactly  but  the  shape  of  the  second  object  is  quite  apparent  in  Fig.  6-16, 
Removing  the  front  object  from  the  diffraction  pattern,  once  it  is 
recognized,  can  be  helpful  for  finding  the  location  of  the  second  object 
in  a  case  like  this  (see  case  j,  below). 
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Figure  6-15 
Amplitude  Contour  Plot  of  Reconstructed  Diffraction 
Pattern  in  Front  Plane  of  Double  Specular  Object 

(Case  i) 
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Figure  6-16 
Amplitude  Contour  Plot  of  Reconstructed  Diffraction 
Pattern  Near  Back  Plane  of  Double  Specular  Object 

(Case  i) 
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Figure  6-17 
Amplitude  Contour  Plot  of  Reconstructed  Diffraction 
Pattern  in  Back.  Plane  of  Double  Specular  Object 

(Case  j) 
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Figure  6-18 
Amplitude  Contour  Plot  of  Diffraction  Pattern  in  Front  Plant  of  Double 
Specular  Object,  Obtained  After  Removing  Back  Object.   (Case  j) 
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Case  j :   A  final  test  using  computer  generated  diffraction  patterns 
was  performed  with  an  input  geometry  as  in  case  i,  with  the  difference 
that  the  distance  between  rectangles  was  set  to  15. OX  and  both  were 
taken  to  be  specularly  emitting  in  phase.   Distances  to  the  observation 
plane  and  from  these  to  the  first  and  last  plane  of  interest  remained 
unchanged  as  did  the  number  of  investigation  steps.   The  computer 
recovered  the  location  only  of  the  back  rectangle  at  115. 25A 
from  the  observation  plane  (see  Fig.  6-17).   This  was  accomplished 
using  7.5  min  of  computer  time.   It  is  surprising  that  the  plane  of  the 
nearer  rectangle  was  not  found.   In  a  case  like  this  it  may  be  helpful 
to  remove  the  clearly  recognizable  object  from  the  output  by  setting 
a  portion  of  the  reconstructed  array  to  zero  and  to  use  the  remaining 
diffraction  pattern  to  try  to  recover  the  rest  of  the  composite 
object.   This  was  done  here:   The  clearly  indicated  back  rectangle  was 
zeroed  out  and  the  remaining  disturbance  propagated  forward  25\.   A  new 
volume  of  interest  was  defined  to  extend  from  5X  in  front  of  the 
new  observation  plane  (as  the  diffraction  pattern  there  was  not 
displayed,  it  was  not  on  actual  observation  plane)  to  40X  ,  which  was 
investigated  in  25  steps.   This  took  6.3  min.  of  computation.   The 
computer  recovered  the  front  object  0.15A  in  front  of  its  actual  location. 
The  resulting  diffraction  pattern  there  is  shown  in  Fig.  6-18.   The 
fact  that  the  center  portion  of  the  rectangle  is  distorted  is  attributable 
to  the  removal  of  the  center  portion  in  the  back  object  plane. 

The  obtained  results  from  this  series  of  tests  imply  that  automatic 
reconstruction  of  objects  from  their  diffraction  pattern  is  feasible. 
Of  course  for  this  to  be  accomplished  not  only  the  objects  but  also 
their  reconstructions  must  meet  the  criterion  which  is  applied  to  make 
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the  object  -  no  object  decision.   This  is  the  reason  why  in  the  above 
tests  certain  objects  could  not  be  recognized  by  the  computer.   Neglecting 
the  evanescent  wave  destroys  a  larger  or  smaller  amount  of  edge 
information,  depending  on  the  object.   This  has  to  be  kept  in  mind 
when  judging  the  merits  of  the  proposed  scheme. 

2.   Testing  the  Reconstruction  Scheme  with  Experimental  Data 

As  the  input  for  the  final  test  a  set  of  data  from  an  actual 
diffraction  pattern  was  used.   This  diffraction  pattern  had  been  observed 
by  C.  Griggs  (see  [Ref.  27])  for  a  different  purpose.   The  object  giving 
rise  to  the  diffraction  pattern  was  a  square  transducer  of  2  inch  side 
length.   The  active  area  was  1.75"  by  1.75".   The  pattern  was  sampled 
with  a  linear  transducer  approximately  73A  away  from  the  transducer 
by  use  of  a  scanning  system.   These  data  suffer  from  a  few  minor  details 

(see  [Ref.  27]),  the  most  important  of  these  being  that  only  one  half  of 

the  diffraction  pattern  was  observed  to  conserve  time.   It  was  therefore 
necessary  to  fill  the  rest  of  the  input  array  by  a  mirror  image  of  the 
actually  observed  data.   This  explains  the  unnatural  symmetry  that  can 
be  observed  in  inputand  output  (Fig  6-19  through  6-24) .   Another 
drawback  of  this  set  of  data,  which  has  to  be  taken  into  account,  is  the 
fact  that  the  difference  between  adjacent  sampling  points  is  0.9\, 
i.e.,  the  sampling  was  performed  at  a  rate  very  much  below  the  Nyquist 
rate;  aliasing  effects  are  therefore  to  be  expected,  mostly  at  high 
spatial  frequencies  which  means  that  distortion  of  edges  by  enhancement 
or  deemphasis  is  probable.   But  the  most  important  reservation  against 
these  data  is  not  due  to  the  manner  in  which  they  were  recorded  but 
rather  to  the  object  which  generated  this  diffraction  pattern:   As 
the  transducer  had  to  be  clamped  at  its  boundary  to  get  a  watertight 
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seal  it  is  quite  possible  that  it  does  not  constitute  a  good  object  for 
the  proposed  reconstruction  scheme.   One  cannot  assume  that  its 
amplitude  dropped  abruptly  from  a  rather  high  value  to  a  very  low  one, 
but  would  think  that  his  transition  would  be  more  gradual. 

Nevertheless  the  test  was  performed  with  these  data  as  they 
were  readily  available. 

Amplitude  and  phase  of  the  diffraction  pattern  after  reading 
in  and  reflecting  about  the  center  are  shown  in  Figs.  6-19  and  6-20. 

The  volume  of  interest  was  defined  to  be  between  60X  and  SOX 
from  the  observed  diffraction  pattern.   This  valume  was  investigated  in 
8  steps.   This  investigation  took  about  4.4  min.  Two  possible  object 
locations  were  found  at  74.66 X  (edge-measure  29.927)  and  at  67.68X 
(edge-measure  26.438)  from  the  observation  plane.   At  first  sight  there 
is  nothing  in  the  amplitude  display  (Fig.  6-21)  from  the  location  with 
the  highest  edge-measure,  to  prefer  it  as  recovered  object  to  that 
of  Fig.  6-22,  which  shows  the  amplitude  distribution  at  67.68X  from  the 
observation  plane  (or  for  that  matter  from  the  input  amplitude 
distribution),  except  for  the  fact  that  it  appears  close  to  the  given 
nominal  distance  of  73X.   In  fact  the  amplitude  pattern  for  such  an 
object  does  not  change  very  much  in  the  volume  of  interest.   The 
plausibility  that  the  correct  location  was  assigned  the  highest 
edge-measure  becomes  quite  high  when  the  phase  distribution  in  Fig. 
6-23  (at74.66X)  and  in  Fig.  6-24  (at  67.68A)  are  studied.   It  can  be 
seen  that  in  Fig.  6-23  all  points  on  the  face  of  the  transducer  have 
the  same  phase  (or  very  nearly  so),  whereas  in  Fig.  6-24  a  distinct 
variation  can  be  observed. 
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Figure  6-19 
Amplitude  Contour  Plot  of  Experimental 
Data;  Input  to  Reconstruction  Routine 
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Figure  6-20 
Phase  Contour  Plot  of  Experimental  Data; 
Input  to  Reconstruction  Routine 
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Figure  6-21 
Amplitude  Contour  Plot  of  Reconstructed 
Diffraction  Pattern  with  Highest  Edgemeasure 


106 


Figure  6-22 
Amplitude  Contour  Plot  of  Second  Obtained  Reconstruction 
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Figure  6-23 
Phase  Contour  Plot  of  Reconstructed  Diffraction' 
Pattern  with  Highest  Edgemeasure 
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Figure  6-24 
Phase  Contour  Plot  of  Second  Obtained  Reconstruction 
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The  result  of  this  test  seems  rather  inconclusive,  but  this  is 
due  more  to  the  nature  of  data  available  rather  than  the  proposed 
reconstruction  scheme. 
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VII.   CONCLUSION 

Computer  reconstruction  of  objects  from  their  diffraction  patterns 
has  certain  advantages  versus  optical  reconstruction,  notably  the 
avoidance  of  longitudinal  distortion  due  to  the  difference  in  acoustical 
generating  and  optical  reconstruction  wavelength.   The  main  drawback 
of  computer  reconstruction  is  the  need  of  knowing  exactly  the  distance 
of  the  object  from  its  diffraction  pattern.   For  an  unknown  object 
distance  one  would  have  to  use  many  small  steps  of  backward  propagation, 
display  the  diffraction  pattern  at  each  location  and  have  an  operator 
make  the  decision  as  to  the  object  location  according  to  his  preconceived 
ideas  of  what  an  object  looks  like.   As  the  manipulation  of  data  into 
form  useable  for  output  takes  up  considerable  amounts  of  computer  time 
it  will  be  helpful  if  only  a  small  number  of  possible  objects  is 
displayed.   Automatic  object  recognition  by  computer  can  be  an  aid  here. 

Automatic  object  detection  by  a  computer  routine  can  only  succeed 
for  objects  belonging  to  that  class  for  which  the  routine  was  written, 
since  it  is  impossible  to  find  a  criterion  which  will  reconstruct  all 
objects.   A  particularly  important  class  of  objects  is  that  which 
includes  objects  with  sharp  boundaries  and  edges,  i.e.,  rapid  changes 
in  amplitude  separating  regions  of  high  emission  from  those  with  much 
lower  amplitude.   Therefore  the  reconstruction  scheme  presented  in  this 
report  was  tailored  to  detect  edges  and  boundaries,  without  using 
great  amounts  of  computer  time  and  core  memory.   A  technique  as  simple 
as  the  one  given  here  can  be  improved  in  various  ways;  in  doing  this  it 
should  always  be  remembered  that  any  improvement  obtained  at  the  cost 
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of  much  increased  computation  time  will  reduce  the  usefulness  of  the 
whole  automatic  reconstruction  idea.   The  present  edge  detection  scheme 
could  readily  be  adopted  for  other  classes  of  objects,  e.g.,  those 
characterized  by  a  sharp  phase  change  or  by  a  difference  in  texture. 
Once  different  values  have  been  assigned  to  different  types  of  texture 
it  is  possible  to  operate  with  the  same  edge  detection  routine  on  the 
thus  interpreted  array. 

In  this  investigation  the  spatial  frequency  approach  to  diffraction 
has  been  preferred  to  programming  the  Fresnel  equation,  mainly  because 
repeated  small  propagation  steps  are  possible  when  the  transfer  function 
is  used,  but  also  because  it  circumvents  the  Fresnel  approximations  and 
therefore  gives  better  results  in  its  region  of  validity.   It  was  found 
that  the  necessary  neglect  of  evanescent  waves  can  give  rise  to  very 
strong  deterioration  of  the  output,  most  of  all  for  diffusely  emitting 
objects.   It  is  tempting  to  use  a  lower  sampling  rate  than  the  Nyquist 
rate  because  this  avoids  the  necessity  of  neglecting  evanescent  waves: 
perfect  reconstructions  can  be  obtained  this  way,  but  only  for  computer 
generated  diffraction  patterns,  or  those  obtained  in  an  actual  data 
acquisition  system  from  periodic  objects.   For  all  other  diffraction 
patterns  the  aliasing  introduced  by  undersampling  could  be  expected  to 
be  quite  severe. 

For  cases  where  the  neglect  of  evanescent  waves  destroys  edge 
information  to  a  high  degree,  the  reconstruction  procedure  developed 
in  this  investigation  will  not  succeed  since  the  resulting  diffraction 
pattern  at  the  object  location  does  not  belong  to  that  class  of 
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diffraction  patterns  which  the  computer  accepts  as  objects  because  of 
the  deterioration  of  the  edge  features.   However,  the  test  results 
presented  are  quite  encouraging  as  they  show  that  for  many  situations 
the  computer  can  indeed  make  correct  decisions  as  to  the  object  location. 
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APPENDIX 


Main  Loop  of  Subroutine  ADERIV 

This  computes  the  edgemeasure  of  the  array  AMAG 


Setting  thresholds 
THRESH  and  THRESL 

ASMPL  =0.0 
ASMMN  =0.0 


Horizontal  positive  step 
at  point  under  investiga 
tion? 


Horizontal  negative  step 
at  point  under 
investigation? 


Vertical  positive  step 
at  point  under 
investigation? 


Vertical  negative  step 
at  point  under 
investigation? 


Are  ASMPL  and  ASMMN 
both  zero? 


f 


ASMPL  = 
size  of  step 


ASMMN  = 
size  of  step 


ASMPL  = 
ASMPL  + 
size  of  step 


ASMMN  = 
ASMMN  + 
size  of  step 


to  end  of  main  loop 


(continued  on  next  page) 


114 


(from  previous  page) 


>0.0 


<0.0 


Is  a  positive 
step  adjacent? 


Add  ASMPL  to  highest 
positive  adjacent 
value;  store  result 
at  PHA  (I, J)  and  add 
result  to  STPSUM 


V        ¥   " 


1 


Store  ASMPL  at  PHA 
(I, J)  and  add  to 
STPSUM 


to  end  of  main  loop 


Is  a  negative 
step  adjacent? 


1 


Add  ASMMN  to  most 
negative  adjacent 
value;  store  result 
at  PHA  (I, J)  and 
add  absolute  value 
of  result  to  STPSUM 


Store  ASMMN  at 
PHA  (I, J)  and  add 
its  absolute  value 
to  STPSUM 


to  end  of  main  loop 


(continued  on  next  page) 
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(from  previous  page) 


i 


Is  maximum  absolute  value 
adjacent  due  to  a  positive  step? 


L 


Add  ASMPL  to  highest 
positive  adjacent 
value;  store  result  at 
PHA  (I,J)  and  add 
result  to  STPSUM 


Add  ASMMN  to  most 
negative  adjacent 
value;  store  result  at 
PHA  (I, J);  add  absolute 
value  to  STPSUM 


1 


End  of  main  loop 
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LISTING  OF   COMPUTER  PROGRAM 
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